{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extracting BP values\n",
    "In this tutorial we will extract reference blood pressure (BP) values from an arterial blood pressure (ABP) signal.\n",
    "\n",
    "The **objectives** are:\n",
    "- To extract an ABP signal from a segment\n",
    "- To identify beats in the ABP signal\n",
    "- To identify pulse onsets and peaks in the ABP signal\n",
    "- To derive reference systolic, diastolic and mean BP values from the ABP signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\"> <b>Context:</b>\n",
    "    There are two potential approaches to obtain reference BP values:\n",
    "    <ol>\n",
    "      <li>Obtaining BPs directly from numerics data:</li> The MIMIC Waveform Database contains 'numerics' data, which is the numerical values which the patient monitor derives from the waveforms and displays on the screen. This approach has the advantage that no signal processing is required, but the disadvantages that it involves reading additional numerics files, and that we don't know how these values are derived by the monitor (since the algorithms are proprietary).\n",
    "      <li>Obtaining BPs from the ABP signal:</li> This has the advantage that no additional files are required, since all signals are stored in the same file. However, it does involve signal processing to obtain the values from the ABP signal.\n",
    "    </ol>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_These steps have been covered in previous tutorials, so we'll just re-use the code here._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: wfdb==4.0.0 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (4.0.0)\n",
      "Requirement already satisfied: numpy<2.0.0,>=1.10.1 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from wfdb==4.0.0) (1.20.1)\n",
      "Requirement already satisfied: scipy<2.0.0,>=1.0.0 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from wfdb==4.0.0) (1.6.2)\n",
      "Requirement already satisfied: SoundFile<0.12.0,>=0.10.0 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from wfdb==4.0.0) (0.10.3.post1)\n",
      "Requirement already satisfied: pandas<2.0.0,>=1.0.0 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from wfdb==4.0.0) (1.2.4)\n",
      "Requirement already satisfied: matplotlib<4.0.0,>=3.2.2 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from wfdb==4.0.0) (3.3.4)\n",
      "Requirement already satisfied: requests<3.0.0,>=2.8.1 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from wfdb==4.0.0) (2.25.1)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from matplotlib<4.0.0,>=3.2.2->wfdb==4.0.0) (2.8.1)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from matplotlib<4.0.0,>=3.2.2->wfdb==4.0.0) (1.3.1)\n",
      "Requirement already satisfied: cycler>=0.10 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from matplotlib<4.0.0,>=3.2.2->wfdb==4.0.0) (0.10.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from matplotlib<4.0.0,>=3.2.2->wfdb==4.0.0) (2.4.7)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from matplotlib<4.0.0,>=3.2.2->wfdb==4.0.0) (8.2.0)\n",
      "Requirement already satisfied: six in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from cycler>=0.10->matplotlib<4.0.0,>=3.2.2->wfdb==4.0.0) (1.15.0)\n",
      "Requirement already satisfied: pytz>=2017.3 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from pandas<2.0.0,>=1.0.0->wfdb==4.0.0) (2021.1)\n",
      "Requirement already satisfied: urllib3<1.27,>=1.21.1 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from requests<3.0.0,>=2.8.1->wfdb==4.0.0) (1.26.4)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from requests<3.0.0,>=2.8.1->wfdb==4.0.0) (2022.6.15)\n",
      "Requirement already satisfied: chardet<5,>=3.0.2 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from requests<3.0.0,>=2.8.1->wfdb==4.0.0) (4.0.0)\n",
      "Requirement already satisfied: idna<3,>=2.5 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from requests<3.0.0,>=2.8.1->wfdb==4.0.0) (2.10)\n",
      "Requirement already satisfied: cffi>=1.0 in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from SoundFile<0.12.0,>=0.10.0->wfdb==4.0.0) (1.14.5)\n",
      "Requirement already satisfied: pycparser in /Users/petercharlton/anaconda3/lib/python3.8/site-packages (from cffi>=1.0->SoundFile<0.12.0,>=0.10.0->wfdb==4.0.0) (2.20)\n"
     ]
    }
   ],
   "source": [
    "# Packages\n",
    "import sys\n",
    "from pathlib import Path\n",
    "!pip install wfdb==4.0.0\n",
    "import wfdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MIMIC info\n",
    "database_name = 'mimic4wdb/0.1.0' # The name of the MIMIC IV Waveform Database on Physionet\n",
    "\n",
    "# Segment for analysis\n",
    "segment_names = ['83404654_0005', '82924339_0007', '84248019_0005', '82439920_0004', '82800131_0002', '84304393_0001', '89464742_0001', '88958796_0004', '88995377_0001', '85230771_0004', '86643930_0004', '81250824_0005', '87706224_0003', '83058614_0005', '82803505_0017', '88574629_0001', '87867111_0012', '84560969_0001', '87562386_0001', '88685937_0001', '86120311_0001', '89866183_0014', '89068160_0002', '86380383_0001', '85078610_0008', '87702634_0007', '84686667_0002', '84802706_0002', '81811182_0004', '84421559_0005', '88221516_0007', '80057524_0005', '84209926_0018', '83959636_0010', '89989722_0016', '89225487_0007', '84391267_0001', '80889556_0002', '85250558_0011', '84567505_0005', '85814172_0007', '88884866_0005', '80497954_0012', '80666640_0014', '84939605_0004', '82141753_0018', '86874920_0014', '84505262_0010', '86288257_0001', '89699401_0001', '88537698_0013', '83958172_0001']\n",
    "segment_dirs = ['mimic4wdb/0.1.0/waves/p100/p10020306/83404654', 'mimic4wdb/0.1.0/waves/p101/p10126957/82924339', 'mimic4wdb/0.1.0/waves/p102/p10209410/84248019', 'mimic4wdb/0.1.0/waves/p109/p10952189/82439920', 'mimic4wdb/0.1.0/waves/p111/p11109975/82800131', 'mimic4wdb/0.1.0/waves/p113/p11392990/84304393', 'mimic4wdb/0.1.0/waves/p121/p12168037/89464742', 'mimic4wdb/0.1.0/waves/p121/p12173569/88958796', 'mimic4wdb/0.1.0/waves/p121/p12188288/88995377', 'mimic4wdb/0.1.0/waves/p128/p12872596/85230771', 'mimic4wdb/0.1.0/waves/p129/p12933208/86643930', 'mimic4wdb/0.1.0/waves/p130/p13016481/81250824', 'mimic4wdb/0.1.0/waves/p132/p13240081/87706224', 'mimic4wdb/0.1.0/waves/p136/p13624686/83058614', 'mimic4wdb/0.1.0/waves/p137/p13791821/82803505', 'mimic4wdb/0.1.0/waves/p141/p14191565/88574629', 'mimic4wdb/0.1.0/waves/p142/p14285792/87867111', 'mimic4wdb/0.1.0/waves/p143/p14356077/84560969', 'mimic4wdb/0.1.0/waves/p143/p14363499/87562386', 'mimic4wdb/0.1.0/waves/p146/p14695840/88685937', 'mimic4wdb/0.1.0/waves/p149/p14931547/86120311', 'mimic4wdb/0.1.0/waves/p151/p15174162/89866183', 'mimic4wdb/0.1.0/waves/p153/p15312343/89068160', 'mimic4wdb/0.1.0/waves/p153/p15342703/86380383', 'mimic4wdb/0.1.0/waves/p155/p15552902/85078610', 'mimic4wdb/0.1.0/waves/p156/p15649186/87702634', 'mimic4wdb/0.1.0/waves/p158/p15857793/84686667', 'mimic4wdb/0.1.0/waves/p158/p15865327/84802706', 'mimic4wdb/0.1.0/waves/p158/p15896656/81811182', 'mimic4wdb/0.1.0/waves/p159/p15920699/84421559', 'mimic4wdb/0.1.0/waves/p160/p16034243/88221516', 'mimic4wdb/0.1.0/waves/p165/p16566444/80057524', 'mimic4wdb/0.1.0/waves/p166/p16644640/84209926', 'mimic4wdb/0.1.0/waves/p167/p16709726/83959636', 'mimic4wdb/0.1.0/waves/p167/p16715341/89989722', 'mimic4wdb/0.1.0/waves/p168/p16818396/89225487', 'mimic4wdb/0.1.0/waves/p170/p17032851/84391267', 'mimic4wdb/0.1.0/waves/p172/p17229504/80889556', 'mimic4wdb/0.1.0/waves/p173/p17301721/85250558', 'mimic4wdb/0.1.0/waves/p173/p17325001/84567505', 'mimic4wdb/0.1.0/waves/p174/p17490822/85814172', 'mimic4wdb/0.1.0/waves/p177/p17738824/88884866', 'mimic4wdb/0.1.0/waves/p177/p17744715/80497954', 'mimic4wdb/0.1.0/waves/p179/p17957832/80666640', 'mimic4wdb/0.1.0/waves/p180/p18080257/84939605', 'mimic4wdb/0.1.0/waves/p181/p18109577/82141753', 'mimic4wdb/0.1.0/waves/p183/p18324626/86874920', 'mimic4wdb/0.1.0/waves/p187/p18742074/84505262', 'mimic4wdb/0.1.0/waves/p188/p18824975/86288257', 'mimic4wdb/0.1.0/waves/p191/p19126489/89699401', 'mimic4wdb/0.1.0/waves/p193/p19313794/88537698', 'mimic4wdb/0.1.0/waves/p196/p19619764/83958172']\n",
    "\n",
    "rel_segment_no = 8 # 3 and 8 are helpful\n",
    "rel_segment_name = segment_names[rel_segment_no]\n",
    "rel_segment_dir = segment_dirs[rel_segment_no]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Extract one minute of ABP signal from this segment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_These steps have been covered in previous tutorials, so we'll just re-use the code here._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Metadata loaded from segment: 88995377_0001\n",
      "20 seconds of data extracted from: 88995377_0001\n",
      "Extracted the ABP signal from column 3 of the matrix of waveform data at 62.5 Hz.\n"
     ]
    }
   ],
   "source": [
    "start_seconds = 100 # time since the start of the segment at which to begin extracting data\n",
    "no_seconds_to_load = 20\n",
    "segment_metadata = wfdb.rdheader(record_name=rel_segment_name, pn_dir=rel_segment_dir) \n",
    "print(\"Metadata loaded from segment: {}\".format(rel_segment_name))\n",
    "fs = round(segment_metadata.fs)\n",
    "sampfrom = fs*start_seconds\n",
    "sampto = fs*(start_seconds+no_seconds_to_load)\n",
    "segment_data = wfdb.rdrecord(record_name=rel_segment_name, sampfrom=sampfrom, sampto=sampto, pn_dir=rel_segment_dir) \n",
    "print(\"{} seconds of data extracted from: {}\".format(no_seconds_to_load, rel_segment_name))\n",
    "abp_col = []\n",
    "ppg_col = []\n",
    "for sig_no in range(0,len(segment_data.sig_name)):\n",
    "    if \"ABP\" in segment_data.sig_name[sig_no]:\n",
    "        abp_col = sig_no\n",
    "abp = segment_data.p_signal[:,abp_col]\n",
    "fs = segment_data.fs\n",
    "print(\"Extracted the ABP signal from column {} of the matrix of waveform data at {:.1f} Hz.\".format(abp_col, fs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Derive BP values\n",
    "Derive BP values from the ABP signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detect beats in the BP signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Import the functions required to detect beats by running the cell containing the required functions below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Detect beats in the ABP signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Detected 30 beats in the ABP signal using the d2max algorithm\n"
     ]
    }
   ],
   "source": [
    "temp_fs = 125\n",
    "alg = 'd2max'\n",
    "pks = pulse_detect(abp,temp_fs,5,alg)\n",
    "print(\"Detected {} beats in the ABP signal using the {} algorithm\".format(len(pks), alg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Make a plot showing the ABP signal and the beats detected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABfMElEQVR4nO29e5hdRZ0u/Na+dXc66ZB750ouQIBAkCRAggpBjh+ghxFHZw6Ojugo6ifHwdtIeLjMjHJzBs84+hnFeWRAHrwdHZVxxAGDEDQJEMAkxFyBkFuTdMi1k77t3vX90bt2qlevqvrVWr+1d/dmvc+TJ53ev9Sq61tvvVW1tpBSIkWKFClS1Bcytc5AihQpUqTgR0ruKVKkSFGHSMk9RYoUKeoQKbmnSJEiRR0iJfcUKVKkqEPkap0BABg/frycOXNmrbORIkWKFMMKzz///AEp5YSwz4YEuc+cORNr166tdTZSpEiRYlhBCPGa6bPUlkmRIkWKOkRK7ilSpEhRh0jJPUWKFCnqECm5p0iRIkUdIiX3FClSpKhDpOSeIkWKFHWIlNxTpEiRog6RkvswxapVq7Bhw4ZaZyMFM44ePYr9+/fXOhsp6gB1Re4rV67EkSNHap2NWLjnnnvw3ve+1xn31re+FfPnz69CjlJwoFQq4fjx4864uXPnYtKkSVXIUXK44447IISodTbe9BBD4cs6Fi1aJOPeUO3r60Mul8OcOXOwfft2ppxVH2pQuNqFGpdiaODGG2/EN77xDfT09CCfzxvj6qFdVRmKxSKy2WyNc1PfEEI8L6VcFPZZ3Sj3UqkEAHj55ZfdwatXA3ff3f93ihRVwIMPPggAOHHiRI1zkjxyuf63mlBWKimSQ92QO1nprF6N3ksvhbz1VuDyy6tG8Bs3bsQtt9wyrBVZiuhQar27u7vGOUkeI0aMAPDmIPdDhw5h/fr1tc5GKN585P7kkxC9vRClEtDTAzz5pDH0xz/+MY4dO2ZNrlQqYeXKlc7HXnnllbjrrrtw4MABWj5TVLB//35s3Lix1tmIhUKhAADo7OyscU6SB5XcP/zhD+Piiy+uRpYSw7Rp03DeeefVOhuhePOR+9Kl6AHQCwCFArB0aWjY+vXrce211+KTn/ykNblvf/vbuPTSS/HLX/6SlD+qcuvr6yPFuZDJZHDXXXexpMWNw4cPY8+ePc642bNn45xzzqlCjpKDUu5Ucu/t7bV+3tPT44zhhpSSZCtRyf2hhx7C6mFujVJttl/+8pfYvXt3wrkZiDcfuS9ZgssB3A4AK1YAS5aEhqlB6PLw1efbtm2zxjU0NACgk7srbjGAZYDTVpJS4pZbbiE9s9r41Kc+hWnTpjkJj7q8//jHP47Pfe5zHFkj48iRI7j55pvR09NjjaOSu2rXnqeessY1NDSQTkt1dHQ4hYKUEvfdd5+TqB566CE0Nzdjy5Yt1rgl6C+DeOYZZ/6GKjo7O7FixQpyvGuiveaaa6qv8KWUNf+zcOFCGRfHjx+XAGR/kcwoFoukuLVr10oAcsGCBda4m2++WQKQd9xxhzXurLPOkgDkhg0brHGLAbkMkEd+8xtz0KpV8jggewFZamqSctUqYyilrKVSSX7zm9+Ux44ds8ZxQ+Vt9+7dpDhqetXEjTfeKAHI733ve9a4D8ycKZcB8tl//VdzkN6ujY0s7QpAfvSjH7XGPf744xKAvP766+1l+MAHJAD50EMPWcvQmcnIXkAWGxpil6FW+NjHPiYByI0bN1rjVBmOHDlCiuMGgLXSwKt1o9ypoC5l1Y6/K56qyC/s68MyABmbmlm9GisAfAXAyGuuMavyJ59EAf3ftFLq7LTuG1Dw5JNP4jOf+QxuuOGGWOlEwWIAha99jbSxLau4GS2lxDe+8Q0cOnTIGqdUsXVvZvVq3P/aa/gKgHk33khqV+nYD6Ks3FR9/fu//7ulBP2risUAFj72mDW9Bd3dWAag8cUXzYk9+STypRJyAERvL6kMxaeftubvC1/4grMMvb29GDduHH7wgx9Y44rFIu699150dXVZ4zZv3gwAzj0yVYZey75bNfvtoAfX+g+Hcu/o6CDNjkePHiXFrVu3TgKQZ511ljXu7rvvlgDkl770JXMQVc3cdZfsBaQEZFEIKe+6y5jecUD2APIEYEyvVCpVVgI2BfXUU09JAPJtb3ubuQxSyjVr1kgAcufOndY4KaXs6elxxiwB5HFAFgEpLSsQVYYXly+3pkdpVymlvOOOO+QqS31IKeWqVaskAPm+973PGvfFL35RApBf/epXzUF33dVfxnKbUdr1uK3NNIXfZ1H4xWKR1P5bHnigkp6xHVatkt25nOwFZE8uZ81bZyYjewDZnc2SysCx+jxw4IAEIE855RRr3P333y8ByFtuucUad9lll0kA8re//a05SCtDTz5vLENfX1+q3ONAEmdHqnJXiqxYLFrjZra1YRmAyTt2mIOoakbb7O2W0rjZq+8bvO+UU4z7BnLVqspKwHbsc/To0QD6NzhteOihhwAAP/vZz6xxX/nKV1AoFJzqaCmAAoAsYD65pK1mzvj0p1mOrt56663OUxrq3kRbW5s1rrGxEYDDS1+6FH3ZLHpR3si3tOufjRiB2wFcXv53KPSVW1eXsT9R23/k2rWV9Izt8OSTyPX1IQcg09dn7sNLluDLl1yC2wFc2tdHKoPr1BoF6rKUa09D7Y20t7db49QeiXUTVC+DZVxLKcl7ZJyoS3J/4YUXjHE+5L4YwMfb280Nsno13vftb+MrAD79s5+Z45YuRbE8uEu5nHVwf6i11T24AeycMgX3AGgwpQWQB5C6Uegi99mzZ2MxgDk/+Ym1k/7kJz8BALxoW74DWJnJuE8uaWXIl/9tAucAUkcXXXbbnP37sQzAqq99zRy0ZAnuveoq3A7gylwOcvFiY2jD0qW4B0DfBReY09NEQK8Qxv4knnqK1P7HFi6spCdN7bB0KUq5HHoB9GUy5j4MYNe0abgHwBpzCehChgglxlztRW3Xczs6sAzAdz7yEXPQ0qXoFaJ/XFvqRJ9kS5ddVjWCd5K7EOJ+IcR+IcRLgd9/RgixRQixUQjxT9rvbxZCbC9/dkUSmQ6DTu6/+tWvjHE6uR89etQY1/THP2IFgM8fPmxWPU8+iWyx2K/IHWrm/7zrXbgdwPWzZllJu/2003APAGGJAU56/bbJqnTJJaRjn6ruXOQ+9+BBrABw1erVViU4ceJELAYw9r77rB255YorcDmAr48ZYz65pBNZ+d+h0BS+LW/UFZ5SbtZTMKtX44P//u/4CoCfHztmLevWceNwD4Cni0XnahCAfdWjrdxuuegiY3/S27+UzxvrruPccyvptT30UHh6S5bgv7/0JbciBzBq1Chz3kPKcFmpZE2PMmmXSiVSnFppqRvDoVi9GneuWYOvAFhhS2/JEnx0+nTcDuCLCxaYy6BNsq69FE5QlPsDAK7UfyGEuAzAewDMl1LOA3Bv+fdnA7gWwLzy/1kuhKjKyyX0QWs7OqcPLNvlo2bKUlVTM1biAfD6rFm4B8ADjmNkSoG4iFaV16ZA5OLFlQG09dvfNts35bQ6Ojqsz5yydSupky7o7sYKAKd9//vOW8BrAHyrpcU8MJYswf8aNw63A/jCeeeRlvm2vOnHAm0kq1Zu796wAXLVKuMz1eSeB9D3xBPW9BSsbVZuC5el1T5nDu5BvxgwoXTRRZX2//H111vbfw2AewDsmTHDmN7+8jPXwH4PQ39pmK2OnxHCrfC1SVta+lLmmWdIk3tDQ0NlEug0tdeTTyInJWm1+GJjI+4B8EJ5RRCG0tvfXplk+7LZ2KsUKpzkLqVcCeBg4Nf/L4B7pJTd5Rj1jtL3APiRlLJbSvkqgO0ALmTMry2flZ9t5K4rXVvcofnz3ap3yRL81+c/X7FRehYuNKanPNzZs2cbY4CTg8Y1uCnkXiqVKoP2Py27/roneOQ3vzHGvX7mmSfrxKIEFx471u+lS2m1A6hEtm7ECNwD4LnyCaZQaAq/aFkiF4vFSlkfvf12Y3KNL77oJpWlS9GnTe7HF4W+v6nyXAUKubtsAxVn85j19t9n6Xf62LF50aoPA8DevXtJ6dnadmm5jRYsWGCMoe4v5P/wB9LkPvHllyvt2vjud5vblbJHAlo79F14YWWS/d4HPmBdpXAiqud+BoC3CyGeEUI8JYRQBuFUALu0uN3l3w2CEOITQoi1Qoi1rs0NCqKQu+3SxqEzz6w0yOsPP2xskLaZMyvqw3YcTg0M10URqneo0rPZBnqd2N5EOGLdukqHH/Xe9xpVz/45cyp18tr99xvr5JUZM9weLmgDQ4+z1t2SJXhXoYDbATx03XVmq+IPf6iU9ap77zWWVV+5GTfBlyzBLz/zmcrkfnDuXGP2fJU7ldyp7W+zIPU4Wx1Tx5geRymr1foi7i90XnTRSQvKoo7Hrl/v3gRdsgT/3zXXVNq1dNFFzjLYJjF9ZbRt/HhjHDeiknsOwBj022F/B+Anon8tFvYS51CTU0r5XSnlIinlogkTJkTMxoD0Kj9zkHtfXx9J9eiD1kbupPPQOEna1MFNGTwAcPBgcPF1EvppCeFQ2qpODp91ljE9fRLoffRRpx1EXaW4JsZnMhncg35/2wR9k9F26kNfudk2wfeeemplcrcRqN5PbOVNqv2541xk5hNnLavmzdtUb+db3lKJ+/0//qMxbt9ZZ50UHpbVp9ojWQO7+KCuoBRcfZ0TUcl9N4D/KB+1fBZACcD48u+na3HTAJjXb4ygkrs+yFzkrmD7AhC94SjK/fjx41a/kqrcOQf30QULSESml9WWXqFQqEwCXeef78xfT0+PtU5UnOukk/rc1q7dS5aQ/M+Dc+dWyOL5r37VSBZURZ6ULUNV7tykzTlZuF7d8MrEibgHwOYxY4wxugW1Y/JkY1z7aadV2nXjv/6rsV31PFFWKZQJAKjuW0GjkvsvALwDAIQQZ6D/uPIBAI8AuFYI0SCEmAXgdADPMuTTC67OomBrEH0w2sidOrh947jJ3VYn+mmJJ2+7zWxpEMk9KbJwEZmqO9tg1BXeD/7mb4xlLRaLpE1GvZ/Y8ldrW2bYKnec7HeUlTFgH696u9pW5LqQoNQxVblXk9wtO1T9EEL8EP33TcYLIXYD+HsA9wO4v3w8sgfAdeXbUhuFED8B8CcARQA3SCl5Xm/oQJQOyqHc9Tifwa3enBeE6gi9vb0olUrIZMLnX+7BvQb9S9CzTz2VJT3fuM7OTmOdUMpKtdvU4F4D4OaxY61xCtSVIKdyL5VKKBaLlddgmOJqRe6ccRwTGdX6oI5XXyHDMdlxw0nuUsoPGD76kCH+TgB3xslUFFBVKjWOOhijdBbqJNDT01M5lxsEp3KnxkVR7pxKkJo36lFYKhlzKPKodeci97r33IlxUSZZjjERtBbDvlawVsq9Lm+ocjRaFEXOQQKcy3fuOqk1Cdg89yiTNsfkHoUsqmlVUdIaynYbNY5aBup4pfYnHaY6qZVyf1OTu40suMk9qnI3Iallea2UO2VDqq+vz7jxmiS5c3nplKvvvoTMPWlzb6hS0isWiwP6limOg7T1duUQC5Q+nCr3mEjSluEmbU4V1d3dPaBMYTGuvFGfmaTCo8aZBiSVeLiVuw+5NzU1keKo6QFD20bjfi5HWantryPuxDjcTssMOSRpQVRTuZdKpYpvR+nwUkrjFe9a1UlSyh2gkXs17RafjVK1f1Ktzchab5QCtDP91PQ4T6MBfH3Y9f24qS0TE1E6HseyjNtzj6LwTM+Noty4laCrM6sNQ6oSNOWPe0WWBFkkQe7VPC1DPZHC3e+SsmW4hIxqV1OdpLZMTKjGyOfz7MutapM7hQQoHSZJpc0V5xoYCmo1QyF3DkLRwTHhASCV1ZcYh2q7KjVbrY3XWm2o6n2Y0jdT5R4BqgIbGxtZNlR1cNktlDgfhafeF0Mh96Go3CkDQ8W5XnFcS2vJ1Q4qbjgod649Es79BcorGaLsaXH1E9U3XStoIUSq3OOgoaGhqv6y3lms71b3IEYqCfgohlp56a4418Cgxg2HMvi2q89lJ1sMNS2fOOqkHddu0+O4xBM1jnq6yfU9yuqZTU1NKblHgarAQqGAnp4e0gkS7qNQcQlKxVHVrIsEuAkvCeWeBLlX++x/NptFLpeLTQIqzlfhV3vCy2az5LJSx1i1NlRrMV7VuGlqakptmShQjaYaxKVmMplMzUjbNXtTSYCqGChpKVDUUaFQcKanNkq5VK9SURQrjdKu2WyWFNfQ0OAsgxCCVCfZbJb03CTInZIWZYIC+q1P16St2stVVqql5UrL1/qgtAO38GhsbEyVexQEyd1V0UrhU9JzbbwpIuPuLHHVUVK2jKuT+izLfcvqaleqLUeNo+zhCCG8JgEucnf1O1WGfD5PWpG5bAOqvSClRCaTQS6XY53IKMq9qamJrf1994NM+VP129jYiFKpZH0DKifqltxdarahoYG0ZKR0Au5BS41zqVluz9WHBKjLcm51pEjWZctRyJgap8id2k9cdUI9LUXt65TJWMVRPHJKGahjgtr+AJ3cKWWllIHLltHrF6jecci6JfdqKfcoHZmL3LnLaktLj6PUST6fhxCCTR25lvn6AJJSOtURt3J32TIAEhMBLq+XSu4UwgNqI3gA++snfFcflBU516pSbwdbHDfqjtyps2gtlTtXh/chPM6TQZS641ZuVAsqCfuGQtrUlQCnCPApa29vL2k1QxkTPsqdY5WibB6Ar69TBYrtmSrOxy2wxXGj7sidQhZA9ZW7njdu5W4aQFTlppM2NY6L3KmbjD5EBtAGGpdyp6bHrWZ9ScV1yIBzv4JS1lKpRD4P72N9cCh3NaFQLkX6bKja4rjxpiN3qvpQ4By01BMESfjQFOVGJTJOz5Wr7qKs3Dg9d8ppGSGEkyx0wuPyerlWM7VauVEmgSTInasMqS0TE0l57pwEpZQAx1FIfUONalVQ4iiDO5/Ps5M718kgHxXFqVI5T8v4Xp7hInfu1axrIqOuUnzK4DNpc7Ur9b5JuqEaEUl6rtSOTCE8n81IrsFN7XyUjqwmKEpZKYObSx1x2zIKPqdlqOTuqrtMJuM14XEqdy7PXZUhLrn79mHqpjCnaEuVe5VAHdy1OBngs3ynkAXVquBSs+o7XZNQ7pQ9EmoZKHENDQ3WM8c+Cp/bgqK2v09fB6oneEqlEnki496MrIVlqM70UyeolNw9USvlDtCPuPkoPJ/BTbUq4qpZNWh9yJ1rNUO9oZqEmuVW7lRSoZ40oZY1ri2XxJ6L7xlx6l6aaX8J8BuvHHH6XgWQ2jLeiDJopTSfh47S4aut8Lh26X1tmVwuV3V15OMbA3zLfK4bqoAfqVDqmJsYGxoa0NfXZ/zKuyRWs9T28lmlUMa1zwQV9xJbqtxjIsrmESWO25ahKrxcLmc8uqbiqDdUOW0ZpdxdeRsOExk1Pa4bqgBdBLjaX99QpypGapyrP3FbH7YxEbVdbWXwzVtcuzX13GMiinIHeJaqVN/Yp7O4rA+9DJwbqtQTP9U+LcO9v8Bl3wHwUnjUTWZKHfv61Vx14lOGuGX17cMUwVNtW4ZaBm686chdgXLmXKXH1QnUcynLd5s6TpLIXMcvfT13buXOub8A0JW7ycP19ea5lHsUz50jzmcvhWN/IYrnDrjJnXI5iWv1kSr3mIjizQG0JSinSqEu322eq66gKGXw8ZeLxaLVc01CuVNXDLabrL6+JuftTioZ+xyFpFhfSUzurjjfI65xj0wGiZFj0vYVY9Tx7xJjKblHRNQZnrostym3JI5CUjqL67Wq+vvXAff1c1V3ttcZJHHOnXKxi7r05SY8Cqlw7kOo1ZFrQ9WnDL77UDaxoMiYY3LXV4LUMtTClqGM12w2a2x/6gTFjboj96TUjEu5UTs8lQQoyt2leqgdnnpUKylbJpfLOU83UJfvtfChfZU7lxKsdlkpZKyXgdOW8dnsdcUldQCCcvIOSJW7N5JQKdT0uL1U1+DWv3WGMjCUfUMlC5dy4z7nbquTILlX+7QMhUCTIHfK5E5dfXJuRuvt6lrNckwCvsRIKWsS45XSh1PlHhPq22m4NyNtDeKzfKOeEacMbuqSlkruLg+/1srdZ/nO0a7AyTqx5U9vLxPhqTJwKnfqys13M5riuUtpP0vuW1auCYpbuVPJnWLLpMo9IlTjKnKnXGIA+E5fUO0Wjjhfq8Jn8xhwW1Dc5G4bGNSyKvh6s67nUsSCqhMAzgtArv2FJGw5ztWsagdXer79pFq2DHCyLxWLRefqwyY89PSGpS0jhLhfCLFfCPFSyGdfFEJIIcR47Xc3CyG2CyG2CCGu4M6wCVEtCGqHj3taAoB3XLWVu88xslp57q6yUid339UMpf0pcRyTe9CW47oARInTJzJXP8lms1Zi9LFluC0oal+nKHcApHHt6pvcoCj3BwBcGfylEGI6gHcC2Kn97mwA1wKYV/4/y4UQWZacOqCrGSFEbH9ZwWdZTu0ElBeW+Sh3blvGVg7fZXnckwbB51arrL51QiX3bLZ/ONgUvs/k7vNOk7hxPspdHV2lkLuP3cah8H1WHxy2DJWbuOEkdynlSgAHQz76FwBfAqCva94D4EdSym4p5asAtgO4kCOjLugdnrKUos7w1EHLsVEaTI8yuG0DKKotw3G+nnszSuWPi9y5Fb5Kj6JmXc/1seUo7zRJ4pw74N6b8SF37s1jLnJ3CQ+fPkypE05E8tyFEH8GYI+Ucl3go6kAdmn/3l3+XVganxBCrBVCrG1vb4+SjQEIViDVguAgC31CoXh4nJ47hdy5NlR9LQjXJBvFlnGRAHVDvVa2DJXcqZM7JY77jLiPLVMqlZxjwtaH1aoil8shm82yTAJJKHdqH6bYN1zwJnchxAgAtwC4PezjkN+FtqyU8rtSykVSykUTJkzwzUZYeip/pIr2Valcg5v7tAwnuftsqAK0wc21pKUu37PZLIQQVffcqXEUcleX07gnd1ccZTXjO0EBdgvKZd/4tj+XLaOe6bOhWg/KfQ6AWQDWCSF2AJgG4AUhRCv6lfp0LXYagL1xM0kBdXb0XapyD25u5U4pK6ct40PuuVwOUkqSv8y1oUqtk2ord5U3FzFS3ryZVFmp5O67SrFtDHOv3DhtGcC+UUp1CxSGvHKXUm6QUk6UUs6UUs5EP6EvkFK+DuARANcKIRqEELMAnA7gWdYcm/MFwN+WqbZyd3nuqgwcyt339QOUzUM9jstf5jgK6TvQamHLAH57OFwrN06rSiltgE/hc67cfAVKHFvG1y0AMLSUuxDihwBWA5grhNgthPiYKVZKuRHATwD8CcBvANwgpaxKSXwrmku568+0xVFtGZUe9RJLErYMp3IHaEqQW7lXk/D0slZ7Q5VqaVBIytVe+kYppQw+cVRbxjV2uD136sTjMwlUS7nnXAFSyg84Pp8Z+PedAO6Mly1/JGVVcCkyn7yp9FzvguHaXwgqfM4NNSCeVeWr3Gply3B67tzKnYvcfSdtbnLnOhnmu6GqDkoIIQbF6HmrB899yINS0b7KncOHVOlRiSyu1xdVuXNtqPqQQNwz/Qq+REa9eFLt0zKcKzfqZOE7QXGVtZrkrtLyucQEhG8KU1eL+nOHtOc+VMFtyyhwD27bJqPvYPQ9LRF3NaOemYRy5xrcAG1y9z1Bwn2JKY5y11dulD0X6umbpE7LcK3cOO90UG00Wxl8V0bAEPPchwt8fa2kNlSpnSUsjuov+nb4TCaDTCbD9vZILuWukM1mSROej0rlmPCibJRyrWaGoudeS1uGukqp5gQV1XNPyd0TvlZFEpeYKHG29Lg7i6/NU0vPHQgfaEkRWRL3FyhxPnVSLc9dTapUC7Katgx1lZIUudvKGtVzT20ZT/gSXjabRSaTqck5d1NcmHLXl3RhcZy7+VTSTuIoJEAbQEOV3Kmnqnw9d1f7c9yMjOq5D8UN1WpbSypvPp57qtw9kZTqTcJzB9y2DHUjh3NDLakNVc4LQFQSGIrKHYDXDVVTXJTJnbJyq5Utw7nnkhS529pBPZfquafK3RNROnwmk2HdFKLEcVwU8p2gXF+4TLVl9Gfa4pIcQFyeO6WfRNlLiat69VubpvRq7blTbp7qExTlUhTX/gLX6kM906cdUs89IUQd3HE7sv7MuHFhyj1Opwqeh6/22W8OHzopC8qX3GtxicmUXq2tKi5bhvL2SGpZFTgtI8C+cucWFNyoS3LnUgI+y22OUxXBZTQlLglbptqeO/eGKtckwLki47Q0aqncXZYRtQzBOC5bRggBIewvjktiP0jVb9geiUKq3CMiqnKjLsu47Buf0zKm5+p2i8/V7WqflklCucclMgWK2opC2tW4ABSl/TnIPfj6gaF4icl3Rc51FFKlZ9sjA1LPPRaS2qCpxsZbcD8ACO98vq8foA7uJF75ayqDHsfhLytQ6wRwb4Jx2W16ehy2ge/xQKqvndQxwmpuqFJFG/eqksO+4UbdkLvv7KgamPI6WqD6l5h8luVUCyLuu2VUWkmdluGwIAD6a1q5bBn1TMBddxzESK0T30lAiQquCcpnwuvrC/+imyirFEpcEpeYTOmlnntM+KpUADU5LeO7zOPyXH2OwnEc50vC11STsY0EuDx3hVochaQu82tlaSS1Sgkra5QJqha2TKrcEwS1ApPqyFyXnfS8AdXdUHMp/KSUm4/qAeK9xCkJ5e7T/tS6s9lyUfeXbKtU4OQ7aKpty3D3dU5y9xEepjiFVLlHRBTlxknunJeYALvCS4K0VZzrnTZJDm6u1QynLaNIlmujlEP1RiE82yo1qI6rdRQyqbJyTlCctgxgv6vBjbok91ood+5z7tzLcq4LQFybUdS4KL4m5zl36iTASXiubzsaKrZM3LKqZ1b7KCxXH45iy6TKPQKiDu44qlKBy3MNKi1KHIdyo8bpzzTlTaXHpWapS19uz13lDaCdquEi9+Bxw6Houfu2K9Vu4zotVcujkKnnngCiDO5MJmP0IfUYgOcSk+8SlBJHUeS1PmlQC1uGa5UCVJfcfT13bnJPwnO3tSv3KgWwvzM96qoiji2jp5cq9whIwpahppeUl0qNMz0z6ZMGrmOkHApPoVaeO1Abck/ClnFtqFItDQ5bJtg3ucvK0Q4ASBOUXieuvbRUuUdAUrYM4HfZhboZSVluA7QOTz0eyDUwbKpSj+NU7gDd0uK2ZVxqy0cJcpK7z4qMasu50tOVtq3ufBU5Na7ap2Woe2lqTFDGdarcI4JagUmRO/X0BYdK1QcQpVNRlVutPPc4JKDnj1u5c1lVAP1r+3xUL6ctZ1OWwf2AOHFhqpdi38Qldz3GVgafvgmAPA5T5R4B1KWPrzen0qvGspzaWYKvHwBorzN4M3juAG3py0nuUawq6kQWhyyi2nJctgzHpJ2EtQj4XyaL87oQhVS5R0SUZblLzdaK3JOKc5GAUkdxCCosrhabzJQbyhQiU3FcVhW1DD4bqhwWJDUuSrtW+zY2py3D2TeB1HOPhCidgPL6AYBGAhyea9TBTY3jqhP1h1u5U1QqtaxxCU/FALUhd+49F9vJsOFAjD59mHIAgtqHfffIUs89AXD7y0Fy51CzKi1TXFSVSk2Pu064FB73RMa1oa6QRJ1wr9zinIf3ecuoai8fYqzFhmo1VmRRxFiq3COA219WMYD9BWMKtbZl4lwAiapmqW/U5K4Tiv/JPZFxkQVlMvbdKHell1RZuew2n4tCXLaMKgN333Tdm0mVewQk0QkUOH1onyNTpvSixHFuqAK8yr3aREaN8yEBV1ywDGF1oivoattynGX1mci4yqDAXVbqeKW0l+uZ3KhbcudUbhQfmsNzVeAmPM4bqj51QlXuSRAZZ/tTb7JS6riWG+rUwwPVsOXCVikckzY1b644/ZnUvAHx7DFu1CW5q2W5a5OJ2uEpnYWijqmdRT3TFeci0KRuqKr8DVXlzrVK8SEBV1zStlzcOvG1ZeKeIKv1DVVVVupkzKHcU3KPiLAOH0buCr4dnmsSqKXn7lKzAJx+NVdZFWphQVDjfEiAEsfV/sGjq670KCu3JGwZDoGS5OqTc0PVdxKwcRMXnOQuhLhfCLFfCPGS9rt/FkJsFkKsF0L8XAhxivbZzUKI7UKILUKIKxLK9yBw+1++FgQ1PR+lRYmrpueq0nLFcW5GUcuq54/blvHpJ1yXmKheL9fKjVpWlS+f9qe2q2sfgoPcVYxvGeKMV/25lD7MBYpyfwDAlYHfPQ7gHCnlfABbAdwMAEKIswFcC2Be+f8sF0Jk2XJLAJf/xa3cg3njIG0fv5JrglLw2WSOs6RV8LG0MpkMpJQkW26oqtkkVm5x9xfUawVccUnupXAodwUfG5XjnDtgt3m44SR3KeVKAAcDv3tMSql60xoA08o/vwfAj6SU3VLKVwFsB3AhY35t+az8zPW6VG7lrmKoeTPFJeFXJkFknIMbQFVPXww3cufaUKV67pS4pMrgmrQ5XyvBfc6dImS4wOG5/w2AR8s/TwWwS/tsd/l3gyCE+IQQYq0QYm17e3vsTAwXq8L3KGQtzjnXisi4j4dSzhxXu07Un6G8oUqdBGxxSb0VEuBdkQ+lPsyNWOQuhLgFQBHAw+pXIWGhOwdSyu9KKRdJKRdNmDAhTjZUeipPXmdTXQQAuDuyAucmI7c3aytDqVTy8lIBv+OhnCcNgHgvcQq2l+vtjCqOWwQMtcld1Rl1EnB9Jyu3peGzwkuC3Klf6kMd19WwZXJR/6MQ4joA/xPA5fJkze4GMF0LmwZgb/Ts0ZGUhwdUf0NVpUUtAzWu2qsUwFx31VDurjift4d2dXUZ43z3K1x1kslkqupX6156LWw5jrKqmFqUQaVFiRvytowQ4koANwH4MynlCe2jRwBcK4RoEELMAnA6gGfjZ9ONpJdvXBdAOHbfk1yWU67uA/EHUJJlrcXynSMu6l4KNe7NYEFxlyGJFbQpPW44lbsQ4ocAlgIYL4TYDeDv0X86pgHA4+WKWiOl/JSUcqMQ4icA/oR+u+YGKWXypUB0sqjmhmpSg5azrEnEVduC4G5/12pGgXPC457I9PPVKh8+eYsTV60VWRJlqMV+EBec5C6l/EDIr79nib8TwJ1xMhUFvhVYK8KrhUpN6hILZTXjeulaPSn3uBfAfMsa1b5Rp47iliFO+4eVgdKHbWXlHK9AtL7J4c1zoS5vqPr4ZJwdmdJZVP641KytrMFlPtfNU464JPcXqqnck1KzNsLztW+SIMa4q1lqGXzKyt1e6rkUa4nDm+dGXZK7DzFW05ZR8Nlk5IyT0nyxx2eCcsUlZcu4XkQWTG8okbt6JlC99q8lMSZpQcUpg0rLtwzc1uKQuMQ03JCUcquW6tXzRi0DJc5VJ74TVFyFF1Vpu+JU3oChs6GqP5MaV+29GWpZoxyZrdaei3qmz9csVrsPp7ZMBERd5sdRWsG4alsVnBdFqr0sT2Jwq7ihaMvoqAW5u3xtipDhPDKZxOQ+FFYfqXJPAElsqClwq9m4HT7sSx3iqqMoG2rc/jJ1AHF+4TJAI7xaKfdqnzThOjI7VG2ZqEKG6/1IqXKPgKSVWzVsmaG+LFdpueKiLmkpp5tsZQjGcS3za0HuSZ0gqoXqHc5HIW1xvn1OpWUqAzfqktx9NlTjHucbKsu8uAOD6qUqxJ3wgoTnmgRUWrYyqPS4Xr/KNeGFlTXusT+VHsflmSivn+Cc3Dn6uorxyZvPC8a4V5+pLeOB4aTcfTaZqrGhmtT1cxXnKqstvaTqRD2TGkeZ8Cjtaksv6Ruq+jOCcbVcpSaxoT5UJ6jUlomAKLNjtTYPFZLs8EPFggjGUQaGSq8atwCHiy1T7U3mKH3d9VZIW5xvGSiXnWppywzXL+sYFvAlRqC/0aSMd/bbJ05/rqsTqHJUaxKo5cAAeM9+D7XTF0l57lT7ptqXmPSVYFzVO1TO6nPsB6m09HIlibokd59OQInzORkQ50sCgpMMl1XBNbjVM1Wa3LZMnEtMSdgynGSh8maL87HRVJzLS6fE1XJy53wDZnBcxxEy6pkqj66+xPF6a27UJbn7LN9McSotFTcc1WyUwc39VkiqLRO3rAq+qxRTekm0q0Jcbz6JfYhq9/Uk9hco7arSAuJfduJuB27UDbkr+BKZLU6BsixTcbXyZrk811pNZFzeLOBu11qRey3bnzrhcd9zqObhAWq7KsS1Zaj1q0CN40LdkDvV14qyfHNZEPpzuZblKi7OC6aSVm5xXrqWpOfObcvVitw5idFnwiuVzN9T6tooDUuvFhvg1Diu1SdVKKa2TARE6ci+Hd703CQGLTU9buXGrdx9TstQy0C9oUqd3E3p1Zrcfbx5bmI0bW5SBE/UMgznSZtjguJGXZI7t+fK+R4VjjhVrmq/W4aSN2p6UVYp6g+nSqXGUTfeuCZtatxQW80kZbdRiFH30uOUISyO6xIT5TY2F+qS3H0HdzXUrMobR1zSg1tK8/FQBc6BASSj8DmJDKBdLecgPNdt0STaX3+mraxRVqnVWn0oUNuVekOVUgb1J86RSW7UJbm7Ni1qZcsoUOOqdTOu1n6lSo9L9VJPS/koPGpcHMLTT5Co9Lg9d4rdQi0r5yrF1U98j/NS4zhXnz6v0EiVuweibG5wbagOBc+VS7lRiZG7TqiTgOlLyKMs3zniOFduSRCeSo+bGKs1aUc9q0+N41x9xrVvuFGX5O7TQSlx3BuqtfQrXcrNZ2BwnpZxEaMC1+aWSgtI/rKT78qNk/CA6qre4Wa3ca4+ue5qcKGuyT3OhqpKC0hGuVfztAzX4KbkLZiez8DgtmW4PXfOTcZqEF4c1cvluXOOCd/LidVeafmIsdSW8UCcwR11o0xhKKgZzvdy2OIUuP3KahHeUCF37smd6yhk3MtuvqsP3Q7kblfTZqlu8Zli9Diq8KAcD01tmQiIoty5bBmFoey51/LWJpdf6XquAsfKbShYFRRS4TiCF7Ws1BM/1VqlcgkUbhuNavFxo27IXSEJIqvWhqpeBmp6HJ77cPEr1XNdl5iSsmVcl53eTOfcqSd+qvnK51qvtHxuFKe2jAfCZkeu1w9ks/FfDQz4zfDAm+9K/lBTeFz7EPVE7j5n8H3ikrBlkl5pKXCttLlRd+Tuo9x84+JuMilwEp76wzUwXDfofFczcf3KJHxobltGIYkNOg5yp5SB0v7BU1VDadIOpkUpqxJtwfxFsVtsZdXzl9oyERCnEyRtVSSl3Hx8zSSIrFbHQ2uxv0CNi0MC1WjXavd1rkmbeqbfpwymS2xxbBlq3lJbxgM+m0cK1XpNQRRvzpZeFF9zqNsy1DqJe4mpluSehF/tq8gpdgs1Pe79hWqd+HLFxZmg0ktMCYB7cKu0gOodD0sqTi8D92kJ1z5ELTcZa3FaRsrBezNJrtw4VilRL7HVg+ceFlerMnCjrsk97oaqArdyr2YcJ5EFnwnw1slQusSkPzPuO2i421VXs+qPSfX6KnKfdqVMZNV8vfVQbS8Fl6vADSe5CyHuF0LsF0K8pP1urBDicSHEtvLfY7TPbhZCbBdCbBFCXJFUxoOIMrg54tQzAbOajdIJbHE+xEgldyoJKAxVb5aijqh1olArJViLyT3uxFiNFVm1lbvpslMwLp/PW+NceeMGRbk/AODKwO+WAVghpTwdwIryvyGEOBvAtQDmlf/PciFEli23FiSh3H1JQCkGrt13jpMmFNKmLsujKMGhaMsocJOFiRjjLPPj1h13X3elF8VaopaVe/WpxmuQkMNIu7e315o3lZ4rbkhtqEopVwI4GPj1ewA8WP75QQDXaL//kZSyW0r5KoDtAC7kyaoznwCSUe5JWBrVUG5RvFSft0JS0qOWoVAohA4MBT0925VxPW+Ub1iilsEnzkV41ZzwkurrPhOZa6OUWlaufSMF6nilkjslbjhsqE6SUrYBQPnvieXfTwWwS4vbXf7dIAghPiGEWCuEWNve3h4xG6HpGmfkcn69OzKXh6dQi2V5Ehuq1PQoZcjn8+jp6THGKVDSy+fzAKpD7noZhtrknnS7hqVXK7stCVtGwUe5U78CciiTuwki5HeDj1MAkFJ+V0q5SEq5aMKECbEfrDeIGtymBolK7lHTGwrKrVoDIxhHLSt1AFHqpFAoAKC3V9xJQGEoT+7DwXOPY8sE0+IsA6dyH1K2jAH7hBCTAaD89/7y73cDmK7FTQOwN3r26AhTbrZlPmD23FR6eqOZ0ku6w1drQzXJOKpyp6ojqnI3rQSqXSdJEp4rvaFQVo4LQElcYkqC3F3K3SYUuRGV3B8BcF355+sA/FL7/bVCiAYhxCwApwN4Nl4WaaBWIJW0VVp6XFSFF3VwUzuVbWBwD27KhiqlrME4k+cepe64V25JkEVnZyd27doVGufzGlz13GrYcmEb6sH0orzK11aGJC8x+WyoUk7LuDZUM5lMZVXZ3d09KI4blKOQPwSwGsBcIcRuIcTHANwD4J1CiG0A3ln+N6SUGwH8BMCfAPwGwA1SyuTNJQysaEXwcchdVwLctgyV3EeMGIENGzbgiSeeCI1TA8g2CXASGWWDNqysruNhqgzHjx8fpLaj1J0iglqR+7p166xlGDFiBABgxowZA+KCm4z5fB6HDh0aVFdhk0C1jsK67JtgWRsaGqx7KT7tWq0DEGHkfuLECWffNI3D4L2EQqEQWifcoJyW+YCUcrKUMi+lnCal/J6U8g0p5eVSytPLfx/U4u+UUs6RUs6VUj6abPYH5BMATfUGyf2NN94gxSVN7kGoSeXyyy8f8PsgCZRKJfzwhz/E448/bsxbUhuqLtVjUuRBhZfP59He3o65c+da08tms1izZg0efvjhSM8NU27Usv7ud78jxV155ZWDYvS8NTc3D0onLK6pqQnbt2/HNddcMyAu7KTJoUOHQtPjPAobVncHDhywlqGhoQHHjh3Dc889Z42j7M0ogq8FuReLRUybNm1QenpcLpfDs88+i7/9278NTU+hUCgMDeU+XBBskM7OTtx77734/ve/PyguSNrXX389Vq5caYyzefP6M1Vn+da3vmXNWzabxZEjR7BlyxZrXNiADYtTy/s///M/N5bBV5GtWLHCGtfY2AiANri7u7tx9dVXW+NUW+zYscMapybiD33oQ9Y4IQT++Z//OTS9YJ3ccMMNg+o6LO6OO+7A/v37jXGq7oCB9RK2fFfQCS0Yp/rbI488Yi3rvn378NOf/hTLly8fEKeTNlWgNDQ0AHDvV6j2P+ecc3Ds2DFj3lR6F1448ER01NVsNfYXTH0zeKovSNqq/b/5zW8anwmcHBNJo27JXeHuu+8eFBdsNAD4/e9/74xz2Tc6CYQh2KnOPPNMaxkOHgxeLwiPU/lSAymsDL4d/o477sC+ffuMcWpwL1y4EB0dHYPypjq68hh/9atfDairMGsh+FlYWXfu3Dko/2FxJ06cAADcdtttzrICwEMPPTQozbC43bt3k9ILO96r4nQytBHjkSNHBqURFqcmARupNDU1AQBuuummQf0qrF2vvvpqHD161BkHAK+//roxb3qc3vfCVlpPPfUUfvCDH1jL6kPu1D2y3/72t85nKtgm4+DEr8fp/fvEiRNYvnz5oNUMN+qO3ClxYeSuN6Ap7qc//ak1PT2Nzs5OY970hg4rg69yV9AHkilvN954I/bs2UMqw4YNG4xx+rP0SSBoGegTjiLcsDh9gtDj4tZJ8N+msobVnYIeZ6s7RaAABhBjsP0PHz5sjVPpBcnVFKegrwiCedP73Ne//nVjnF4Pa9asMcbpZdUnC5Nyd8Wpf3/wgx+0ltWH3K+77jocP358UGxwBfWFL3zBOvHoaegTbtBa1NtVh76CAk5yw3e/+93QeC7UFbkHOztAI21gIBkH41Qn+Jd/+RcyMdrUjD4IdXURjPvLv/xLY1n1OIWgcg873QAAX/7ylwelp+L0NKjKLUy564pMQd/bsBGZTc26lrPBOgnzccPay1ZWvZ+EEa6KGzNmTGhcsAwf//jHK5/pZBFVuSuMHDnSWIaw/IbF6aQdJCtT++sTrq2vh1lVqt91dXUNymdYetlsNpSwTeMwuE9imrRtNpq+Yre1q95vTXnTEVYOTtQ9uYep5LBBqxNUMD09rq2tzRindxbboNUJVG/gYNytt96KZcuWARjogZoGd9A2CpuggvkEBioLnSCCnU+fBKjkrpfVRgL6yREbuau9gJaWlkF5C0NYnShQyX3ixImV3wcHsJ6evlFqI4EzzjgD//3f/+2M04nWpiwVbLacDls/0dvVRu6mScCmevW6C67cTJN2mAi477778MADDxjzRhVtetzevebrOH/xF39R+TmsvdSY0OtEL4+pHahuQ1TUFbmb7I5gnEIUcg92Aj1u9OjRoekFO6iuJm1Elslk0Nra6oz7t3/7NwDhZBxGsrblu07uwTrRJwEXuYetGMImMhV35513VjZdbYT3jne8A1/60pecx9LCnhksq543W/vrN6htk8CoUaNC48LypiYnmwj4+c9/junTpw/Kn6msNiLTETZBhZF2cOVgqjubQNHP8kdZkQX7iULQWtLT0sdh2CZo2KRtm6AWL16MX/ziFwAGtmtwgvrRj35U+SxYVr0d3vrWtw56ZhKoG3IP+loKun8LmEk7bEkVpnptg1vvVLaOrFsxQVLR44CThGEb3B//+Mdx0003WQlKT9NWJzpB2YhRH4z6c4MdXl9u2+IaGhrw2c9+FoC97oD+Cairq8tqab322muYOHGitU70drXFZbPZyia5jRjPOecc3HfffQDc5K76ii1u1qxZuPXWWwflLxj37LPPQghhzZsOW5w+advIXVf/trx95CMfCX0u1dII9hMF2x6JPkHZyjpp0iRS3gBURJZtMl60aFFlRRFsV31yeuqpp/D2t7/dWGYu1A25Bzuy2gG3Na5LuSvoFoAtvUWLFlXOads6iz4wwuJ0KCVNIbyenp4BaUcZ3Lq1YCO8s846KzQujGRtZQ1Ts644yoQ3Y8YMXHDBBc4J74wzznCWFQBuueUWtLS0WCd3APjoRz8KIJ5y14mA0v4XXHABrrvuOmcZHn300UFpBeMymUylTmyT+9ve9rbKSTRbO1xwwQXYvn27s070zVbbiRQFncDDyqpgK+uECRNwwQUXDIqztZfNljHFBYVnNpvF2LFjU3KnIjg7Xn755fi7v/s762DMZrOYNWsWAPvgnjlzZuX3tvQaGhrw2GOPAbB3lqVLl1Y+cy23FZFRyB0YvEQO6/C2MmQyGZx//vmD0grGjRkzpnL0y1YGRXamOL3NbGX1JTwV51Kzmzdvxlve8hYnMQL9A9eVXj6fR1NTE5ncbct8wK9OXGW48sorsXjxYmcZtmzZgtbW1kHtr2/QZzIZLFu2DPl8nqUP6/cWbKelFII2TbAM6tCAq68ruyUKufu0VzD/o0aNMp6G4kLdkHuYLTNq1Ch0dXVZ1ez27dtx5ZVXWgdGJpPBq6++CsC8fFQIIx4Fld5VV11V2RiMMzD0Dq4UdzA9XzUDAC+88ILT0jCVNZi3efPmhU4C1IFhi6PUnasMQggSMar0KNZHmMIP5s3UXqayUuokrAxBEgwjlbC4kSNHWid3PY6D3L/xjW/gn/7pnwDQ9hdc7XDbbbdhxowZzpUWVTzZlDsljiIUuFE35B7WQU3LfL2iM5kMaXDPnDmTpNxsnUWHusocxYIIG9yKaG3kvn//flxyySUkgmpubnaSRWNjIzKZTOS8mcrKqdwppE0lMhMxugZuWPtnMhk0NzezWFCqDL29vQM2mk2CJ9j+YXGm9qeSe7CfZLNZa1mz2SwmT54MwN6uH/7whwHYFbkCxUZTkyylLwkhQpW7y5ZJlXtMmDoy4K5oyvJdpUdZljc0NJAGI0AbtD5xNltmwoQJJDWj0gtblrtUb1jewiaBsDg10KhLX0qddHR0DCBXCkGZ4iiTu4pzKTxVDu72p9gB1Mk9inIPay8hBKlOKGV98MEHcf3115P6MKWsSty5yF2VwbahCpjJPUx49vb2Jvoagrohd9MgA9wdnntZHozzGYymOCoxumwZn8FNqZNgnGlgmAaQ3ukzmQx73ZVKpQEndrjbP/hMFUchdwqp+K56uNqf05YJe26cdqUoclNc1LyFpRfWhykbquqZwThu1BW5h/mLwOAKdCkolV6S5B5GxmH5891QDQ60sDqhDu4olkZYhw9Lz3TEzVR3ro1XBV+rylTWsLyZSCUITuVObX+fVSplNRNVucchdzUmbOfmgf76PX78+KCLXVHbKyq5h/Xh5ubmQfYNVXhyo27IPWx2tB2t0xHmVwLuDqrSi9KRs9ksmpqaQgeGjqamJgghItkyJsXQ09Mz6MZr1MEdjDORdlQS8PXwg88Eoq/cgoi6oRpHuVM3Xqmr1FGjRqFYLJLaP47nTm3XsEMBYf3JpY7DnklV+NQVGaVdhRCh4z9V7jFArUBTBwXckwAnuavnUi0NTlsmLC7qaYkguccpq8pfFOVOXc2Y8kaZ8FTeXKqXU7nbREBU5Q64J4EklDvVc6co97CyBhF3vAYxevRo570ElT+K5w6kyp0EagXaOjynX0kldy5i9LEggMFH66jKLVjHccndZRuFKXf1TUZRfeggqHcEWlpaUCwWB73kykTuwWdFJUZKnI9yB2j9hJPcTSeIqLZMMC2AftwweCmKy3M3rVLD4kztkJI7AXFOy8TZjAKiKTLAT/VwkbvP4La9W0bBRO5RSBugbVrZTjgEy0AtKyUuOCBtJKCfhDApQcoms3oul3KnXrJR7eoixmp47sG01KsbgiraZMv6jmsqufvEmfKW2jIExOnIcZV7EKZJIAhO5VYoFAa9DjXu4D5x4oRzcAcVfjU8d1N6UdvVl9zVgLQNbkpckLRtZXX51aZ2DbPbKGVVRBv8XoKwMuiTQBxyVysyygoKoCl3ShznaRlTXKrcYyCsI+dyudCr4HEGt483q8cA8ZbllHO43IQnpXQO7rjWkkvh+5BF0uQetD5c5O6KU+2qPo/TT8KUO9UOsJF7cHMzbLLQ+4mrvWxlzWQyaGpqYlfurrJSFfno0aMHnNLxWX2myj0GTG+FDPP6om4y+bw3IsqgVQjrVJSNHAq5xxnctuV7MG+ussY5Ckktay1tGYCm3PVz+HHIXQmZKBOeibSBwe8XcqVnq5O+vj7nJBDWn+Iocmpcd3d35TUlYSujsOf6CA9T/abKnYCwRgPCd+mjLlV9yL1UKpE6smuzDwDGjRs34FuMTMRIGRjcVlXQvnGRsWvCUwNNrY58JoGok7avcqccwdPjbKSt5y/Onot6LpdACZ6+cpXBFUedGCmTdhK2TFjeqOTusmXChKd6/USq3AkII20gvKJNl51cR+aol6Kog5Zqy4wdO9b6FXW29OKQtk+cmshsvjFlwjNZGhTlHoxpamoKfe0Bt3IPgurNB/tdHFtOpcdFeMGVW7XInSJQmpubkclkyLYM5RKjnjeqcrfZMvopHZPwTPrlYXVD7iZbJmyDjqrcg6B6qdRBq95V4VKz48aNw4kTJyIt3+Mq9+BAM02MlIkMGFwnrvSoyj2sXan7EL53BDhtGUpcnOOGwZhTTjkFwOAvqzbFqW8LqpVyN7Wrj69Nub+ix8W1ZVpaWiClHDAxUriJG3VD7rbZ0aXc1S6975I2bkceN24cenp6nOpo3LhxAE4OSJNiDLN5oqpZVYbgt72HrSqAk9+PSq0TG2nrcbZJwEUCQPTVTFjefAa3HmdKj0qgVFuGotxzuRxGjx49aCVo6nMqLu6kTbU0xowZM+i7dk3j2vQF4gpRj4eayD347Vk+ZaVwEzfqitxNtoxL4Sn/i0u5+ZA74B5AikD1uLDOEqZ6gnFKzbomsvHjxw94pilO5S048QTr2IfIAPckEHY8kEruYRexgMGbjCaFH1W5BzFmzBgA7omxpaUFPT09zpUbhdyBwXs4tvY/cOAAqayKaKmTtilu7NixzlUF0E+0nBuqlLxRyxq2wg/rm6lyJ8Jmy7g2N1RcFEsDcHdkUxyV3MPiwspAVT2UzUj1TDW4TXHBicfU4U1kEdXmUe2lPje1K2WZr17THDz2p38JtMqr3p9MeWtqakI2m3XGUcnd1E+CGDNmzKAvejaRO8WWyWQykYRHWBy1XaMq97C4QqGAxsZGNuUedeVmO8mXKncCqJsWcZbvvhtqrrjgoFWISu5K9egePpXcw5bHQojIyt1F7j6qBwjfUA36mpR2pba/KU6vO1NZg56wKU6tZoLkHnxucKI1lXX8+PHOyVil52rXTCaDsWPHOpU7VaBQyX3s2LE4dOiQczNy9OjRpIksygmiuNZS2F5K6rnHgG3pox+ts5EAl+duGrRRFVmY524atMVicYD/SVH4Ycoim81izJgxZOVOJXcqCbhOJATrzkbaUSb3OOQODFRlqgzBlcCoUaOQyWQqbeFa9ehlDXvm+PHj0dHRMcC+CStD2Okr1yRgKmuhUMDIkSNZyV1KWVHlLiHjKgPFqjL1uTBBIYQgnZbR04uzbxAHdUPutqUPMHDJFdbhw5aDwbhCoYBCoeBUZGowqu8OrZbnTo2jKrzx48cPyluYmslmsyTPPZPJOJUgVUVNmDABANDe3g4gPmlT4yiKXMWpMqgbjWF7PaeccgpZBOh1Z2ovPc40JijKXaXnaq9geq6NV1fdhYmFOGWgHKhQ9ljQWjSd0tEnnrAyhJG7aZI9fPjwgFd8cCIWuQshPieE2CiEeEkI8UMhRKMQYqwQ4nEhxLby32O4MmuDrXGBgZsgYRU9ceLEChkD/QPStPSlqJkxY8Y4yZ3qV44YMQKNjY1kcncNjCjLd5sFoS/fTWomm82GLvOjHoWkkrtOnra40aNHk+LCyN0VZ1LuwEBR4bPqCWsvVSeuSWDcuHE4cuQIisViJS6sDGHtb4pzTQLBl73FFTLjxo1DR0fHgBU5ldyD7ZDL5TBu3LjKeDWRezA9qn1jm2RLpVJi6j0yuQshpgL4WwCLpJTnAMgCuBbAMgArpJSnA1hR/nfisM2OwEnCMw3aCRMmDCD3sE4ADJwEbGpm0qRJzrh8Po+WlhYngQKDB5ppkOlldZG7y5unKjdqnejpmUg7n8+jsbHRSQJhKtU0aR84cKCinm2Ttj7hmeJ0S4NKArY4Crn7rMgAGrkDGGAHRW0vlR6lD1M2o4OTtqsMrr5O3XObNGkS9u3b5yxrWLtSVp+2MgT33LgQ15bJAWgSQuQAjACwF8B7ADxY/vxBANfEfAYJtg4K0Ejg8OHDA669m+KCnSDOJEAdGLrHaCqrj+rp6urCiRMnnHGUvEUhdxcJcNkykyZNQqlUGtD+Ye01YcKEAeRui1PPNNktQD8JqA0/WxyF3AuFAkaNGuVty1BJhSICXO3q24dNcZMmTQKAAWPMJmRcZdDbATBP2mFizNSuLlsmeEqHaqNyIzK5Syn3ALgXwE4AbQCOSCkfAzBJStlWjmkDMDHs/wshPiGEWCuEWKsGSxxQO7KNtAEMGLgu+8Y2aH0ITycoU3pBYoxrywC0wd3e3u4c3LrqsQ2MKORuUkejR49GPp8nkTvgJgu9rK64Q4cOoVgsWif3MCVosmVcZ//Vc6nK3aV6qSJgwoQJ6O7uHvASvDgiQJ8YTXFqHOpjh0qMpmfqfU5KaRRjwfaKY1XpNh919cGNOLbMGPSr9FkApgBoFkJ8iPr/pZTflVIuklIuUiosDkxWhY9yBwYqQdOg3b9/P6SUTuVOIbzW1tZKnJosXGThQ+6UOjGlN3HiRNLg1icy25I2jNxNJ1JcqkcIMWhijEvuvb29zo33CRMmQEqJgwcPWkmgtbUVR48exYkTJ5zK3aVmgcG+tqn9hRBk5e4SAaruVH+35e3w4cMoFouxyX3EiBEYOXKks69Tb21PnDgRx44dQ2dnp7XP6crd1a6vv/76gDiK4LFZxkNOuQP4HwBelVK2Syl7AfwHgIsB7BNCTAaA8t/7LWmwwWRVBE9p2AYtgAENbJoEOjs7B7zX2RR38OBB9Pb2WuMmTZo0qLOEkXuwU4Wl1dDQgObmZtKROcBN7pMnTwYAvP76686BcfToUXR1dZFsGX1iNE0WOhmHxejpqbg45B7cjHTFtbe3Wydjve5c7XrgwIEB/SSXy4WW1aW0wzat49gySvDs27fPObkDGLDyoZI7xSLhKEN7e7uVtCdOnFjpwy5y379/P/r6+qxlDY7X4ea57wSwWAgxQvTn/HIAmwA8AuC6csx1AH4ZL4s0UDu8S7nrlosrzrUsBzBgM880uBVR2CaB1tZWdHR0oKOjA319faFpAYMHkK1TuUigtbUVwECCiqPwdHVsi5sxYwZ27txZyVtYfQTLampXVQaKcgfcHr4+CbhIAADa2tqscVOmTIGUEvv27aucXglr26AdQJnw4hKjPjHa2mvKlCmVsrrI/eDBgwMsLdfkzmW3KkIGzBOKinO1q9rDca3I29rarGU45ZRTIIQYeraMlPIZAD8F8AKADeW0vgvgHgDvFEJsA/DO8r8Th40EdNVDJXfbaRkV51LuKs41CajO4poEgP6BZiN3XTHE9dx1cqcoNxcJ6M+1DYwZM2agra0N3d3dVuVOIfeWlhY0NDR4K3fbqRpgoHIPi9OVuy1OEePevXut7R+0ZWx93UXuo0aNQi6XI6te16Stl4HST/T2p1gaYTHNzc0oFAoDxrVpA1yVgSLG9u3bZ514qIJn8uTJlXFjE56nnHLKkFTukFL+vZTyTCnlOVLKv5ZSdksp35BSXi6lPL38dzLTUgBxl+9qg45iywADyT3OJKB3Kgq5t7W1oVgshi7dVZyL3NWrBajK3aXIqKpHt4NsA2PGjBkAgN27dztVqovchRCDrK84ezO6LeOyW4CByt1m37jIXVlfnZ2dzr5OqZPgyRWT0gboyn3v3r3W9tfrjrqHY9tz0Y+l+qy0XXEuRQ64BU9rayt6e3tx6NAhax8OvueHE3GPQg4ZuCrQNWiFEIOOubnI2LWhquIopO3yZoNKMI5yD/NmTRs+uVyOXbnrS1obue/cudNpyxw+fBi9vb3GdgUGn1yh2jJhdazHuSaybDbLptynTZsGgDbhuSZtgGZ95PP5yuUe1+QuhMDevXsr1lKY+NBVtCs9tZKlEqNpTEQRWb52G2XVa+qbd999Nz75yU+GfhYXw57c29vb0d3dba3A1tZW7NmzB4Bd4QdPfbiWeVy2DDBwEnApBhe5qw0628DQL22Z6iSTyVQmC05ydw3uU089FUA/ubtsGaB/mR+X3EeNGoVCoTDgXHdYXD6fx+jRo502WjabRWtrK3bt2uVUjJlMBm1tbdb0pk+fDqCf3F17LnoZTHU3ZcoU7N27FwBtErC1Vz6fx8SJE7Fnzx7rpnDY2AmLU1alq131fQhTnTQ3N6OpqYm0oaryRiV3W70FRZsp7v3vfz8uvfTS0M/iYliT+8qVKzFx4kSsXLnSWoFnnHEGDh486OwsQXIPi2tsbERLS4uTtEePHo1CoeCcBKjL8nHjxiGXy1UGkGlwK2Lcvn271b6ZPn36gE1LWyd1kbs6vuYi7enTpyObzeLll1+2Ln2nT5+OfD6PjRs3OpU7cFJFU8jdtCJTRytdG2XquS6yAIDZs2fj1VdftfaTbDaLSZMmDSBGm3LftWsX+vr6jO06YcIE9Pb24uDBg9Z+MmXKlIrgsbW/qjtbu6r0dOUe9lzVN1955RWnBQW4LUgKuQshKuPa1g76EUxbuzY3N2PChAmVPmxqe7X63L59u7V+k8SwJvcFCxYgl8vhv/7rv6wVeMYZZwAAtm3bZm2Q4AUlV5yNtPVOZevIo0aNQktLi3NwZ7NZzJgxA6+++qp10F5wwQUAgOeee85KAqeeeip27NgBwE3uuufuIlDbwCgUCpg9ezY2b95sXdIWCgUsWLAAq1evdvrLwEm1ZcubImPX5O463QKcXPW44ubMmTOAyGxxW7duJZH7jh07rIQ3b948AMD69eutcVOnTq2sAm39JNiHbf3p5ZdftsaNHDkS06ZNw5YtW6z2zdy5cwEAGzdutOaNQu56GVztNXHixAH2qCnu3HPPxapVq6x9c9asWRg9ejRefPHFlNyjYOTIkbjiiivw6KOPWkn7zDPPBAC8+OKL1sGtiKyvr89JAi7PXcW5yB3oV3gvvfSSM71Zs2Y5yf3MM8/EqFGjsGbNGuvgPvXUU7F///7K5Q6XcreRMTDwcpctbu7cudi0aZMzbsmSJVi7di26u7utAwjoV4K2gdba2oq+vr7KALe1A4WMTzvtNGzZsoXUrnv27KnYH6a4c889Fxs2bKgQXthzm5qaMGfOnEqcqV3PP/98AP193UaMU6ZMQV9fX+UIpilu0qRJaGtrQ29vLwAzuS9cuBBbtmypeOCmsqr2tyn8M888E/l8HuvWrUOxWDSmpTx3KaWT3F22J9A/yW7bts3Z16+++mps2rQJu3btMsYIIXD++edXyN30zCQxrMkdAObPn49XXnkFPT09xoqeM2cOTj31VDzxxBNW0j7zzDPR3d1dWUqbOsucOXMGDG7qJGCKu+KKK/D0009XVIipHGqZb+vI2WwWS5YswerVq62Ddvbs2QCAzZs3W0lgxowZ2L9/f+V1uKa8TZs2Da+++qqTtM866yxs3bq1QhamOrnooovQ1dWFP/3pT9al78iRI/H8889b2/Xcc88FAGecUtquvJ1zzjnYs2dPpb1Mce985zsBAL/+9a+tcfPnz8fRo0exY8cOZDIZY93Nnz+/oshtRDZlypQKMboU/gsvvGCNmzt3Lo4dO1ZZ5Zni1IrxmWeescbNnz8fGzZsQGdnJ4DwOsnn8zj77LOxbt06a99sbW1FT09PZf/DVCdnnHEGNm3aVHmXkilu3rx52LRpk3Ncv+UtbwEAbNq0yUraCxYswLp166zclCSGPbnPnDkTxWIRbW1t1ll0wYIFFXVsahDV4V966SUA5sZdvHgx9u7di/Xr1wMwd5Y5c+Zg8+bN6O7utsZdeumlKBaLWL9+vXVwz5o1C+3t7Th69KgxLaBfWe7cudM6aNUmzmOPPYbe3l5j3HnnnQcpJf74xz8CMJP2+eefj1deeaWi3Exx8+bNQ09PD7Zv326NU3sR7e3txphMJoO3v/3t+MMf/mBVRxdeeCFyuRz+8Ic/WNv/tNNOQ09PT2Uvwqa0AWDdunWVfIRh4cKFyGQy2LZtmzVOpffiiy9a23XevHnYvn07urq6jO0F9E8+GzdutLb/BRdcgHw+j6effpq0EnjuuecAmEl70aJFAIA1a9YAMNfdggUL0NXVha1btyKXyxnb9rzzzqusUkxpnXbaaQD67VYbub/jHe9AV1cXfv/73wMwt8O8efNw4sQJvPLKK9Y4/darjbRPP/10dHV1Yd++fSm5R8HUqVMB9G802WZR5THaBvfZZ58NANiwYQMAc+NeffXVAIDf/va31rjFixejo6MDmzZtAmC3WwDg5Zdftg5uaty0adNw8OBBHD161Oq5zp07F8888wyKxSLy+XxonBrca9euBQBj3MKFCwH0K0HAPoCAkxOo7Zgb0H9s0tau5557LrZs2YLe3l5j3IgRI7BgwQInuSuSVUTmKoPqJ6a2yOfzmDRpktOWOf300wH0b77Z2nXGjBkolUrYtWuXldxnzpzpnNxHjBiBuXPnYvPmzda4+fPnA0Blcrf532PGjKkQoyluzpw5AIAtW7ZYyzp79mzs3bsXnZ2dzgnl6aeftpL7JZdcgmw2ixUrVgCgt6trI3///v1W0lZ7Qq+//rpx3CSJYU/u6sjR8ePHrRWtzkR3d3cbG62lpQXTp093Dtrp06dj9OjR2Lx5szVOTRYbN260xuknZmxERiV3/bSBLW7OnDnYsWOHVblPmzYN+XweL7/8MgAzuau8vfbaawDMpK32P1wDSL1U6cCBA9Z2nTZtGnp7e/HGG29Y627+/PnYsmWLldwXLlyIxsZGPPXUU9a8TZ8+HYVCwanwgJOTlC1OldWlyJWQee2116xx06dPR3t7O44fP26NO/XUU/Haa69ZrY+RI0eitbUVW7duBWAmbaC/37msD3Wkc+fOnaSy7tq1y3rs95xzzsHvfvc7K7m3tLRg1qxZ2LJlCwBzO+iTrC1OvaCtq6uLRO5vvPEGCoWCMS4pDHty1wePi9yB/tnWNhjnzZvnJB6gv5MqIrORANDfQW1xLS0tyOVy1g4KnCRQwDx4gJPLRtfgnjZtGnbv3m1VbkIItLa2Yvfu3QDMg1tNUOp4namsI0eOxNSpUyunkhoaGkLj1Fef2U4tAQNPzLjaa9++fejq6jLGFQqFiidsK4O68KaOV9ragkLuuVyu8r27trR0JWhrV9XXXXFTp051tiuAyiktV5wSFbY4lbdjx46R0mpra7PGzZ8/v3LSyFZ3U6dOddptY8eORSaTqfRNE5+oS4CAWewAJ9vLFZcU3nTkblNuQP/mi1IptrgpU6ZUBoZtN1+dTbfFqavUthhVBtXRXXEKtoGhThvYbBkAlcs4gLmTNjc3o6WlpTLhUTu9SdE0NTWhqakJgL1dVVr79u1zthfQ75O6Vkcu5Qb017G6CUxV7pRJgELu1LRcKwH9nUuuyULd3LXFKVFhy19DQwNGjx5tjQEGfluUS/C89tpr6O3ttcZNnjy5Yo+Z2iuTyQz4uj2KqGhsbHSWATD38yQx7Ml95MiRlQ5nq0Cd8GyNNnHiROeRRKC/cdXJEFNcJpOpXFBxpUcZ3EIIUhyV3MePH+88vwzQyB3oH0BKHdnaQqlyV3pqwjOpe+DkILPZbQC9/adMmVL58gxXP3GdlgFoyl2Po64CXJM2JU6fLCgrAVccRbkDJycBW4xqe8ozS6USjh07FnsFBfTXifqmJUp/svVN/Zmpco8AnfBsFU0d3HqH54hrbW11XgACTnYE13lYTnKndj71EiRX3OTJkyueq43c1cDN5XLW8qo4mzqitkOU9qfWMZVUKKrcFqPsO2DokzvFNuSoN2o7cPcTFWfrm4VCofJl2Sm5R4TqCNQlEofSonYq5UW74iikDdBIYOTIkZWJjoME9OW2rZPqdUJR7q6lKmUAUdtLLwN3+1dLuetChtquVMLjmAT0Oo6r3PXVHTVvHJaWr7izCUrgZFuktkxEUJS7vszjbFxXnK/Cc5E7lQTUc6s5uPU6jmu3ADRyHzFiBEaMGAGgumQcRbnHJXdq3HBR7ra0crlcRfVyr6Co7UrZw3ORtqrjVLlHBIXc1dv8gH5lawK3EvBVjC5bhqLc9edy2DLUXX/qBpJ+EsYGVQYbuQMnB5qtrC0tLZW8UwmPw+bT60RtENvi1LV8V5ytDI2NjZUJbyh77i5Qyspty+gCxdZeKr2uri5jjB6XKveIoJA7gEqHV4ogDLVW7uo2qwlqALkmAQrhRdmg4yB3Fae+BNsEVQZqWW2TthCiMqmoY4e2vAEDB7rpma78+WzkuWL0OFt76c+tJrlTTssAJ/uw2rg0gXKCzPcEkStOb/NRo0YZ41SduMqQKveYoJK7Ihxbo+mNqyaDMOgd3qYsuW0ZdYvO1VkotoxOclQS4CB3dW3cBUUW6r02rvzZ2lXPk43c9faikrut/c8991w0Njbisssus+ZNtauufsOg6tjWNwEaMeqTIbdyt1kaql1d30BEmaAKhUKlHNRJgLLZD9jrWF14uvjii40xQG09d/v0P0ygbrO5votQdRKbctcbQe/UQeidRd80DUJPwzb5XHjhhQDgfHH/u9/9blx11VX47Gc/a41Ty3t1vjsMukqk2jIcnrt6pYF6eZkJ6tKWrb0AmnLX82Qjd52gqO1qsxlaWlqwfft2Zxkuvvhi3HHHHXj/+99vjVPPUjezTVCKUr/4FoROwK5TUJS45uZm3HDDDdY+B9BtGcoEBfTXcUdHB1l42J6vLh4C9gnqoosuwtNPP115ZYUtb4B9HCYG9Y0ttfyzcOFCGQfPPfecBCC//OUvW+M+/elPSwDyhRdesMYBkADks88+a4zp7e2txPX29hrj1q5dW4lzYfXq1fLo0aPOOAruvvtuCUCuWrXKGjd79mwJQD7xxBPGmI6OjkoZ1q1bZ4xbs2ZNJa6zs9P63A0bNshdu3ZZY4rFovz6178ut27dao373Oc+JwHIm266yRq3bNkyCUA+//zz1rjly5fL++67zxpz6NChSll7enqssZxYt26dPP300+WGDRusccuXL5cjR46Ue/bsscb9wz/8g5w1a5bct2+fNe7222+Xl19+ubWvU7Fnzx6ZyWTkzTffbI37/Oc/LwHIj3zkI9a4MWPGSADy4YcfNsaUSqVKe+3cudMY19bWJgHIKVOm2AtBxBtvvCH/8R//UR45coQlvSAArJUGXq05sUsGcpdSyhdffFEeO3bMGtPV1SWfe+45Z1pf/OIX5dSpU+Xhw4etcf/5n/8pf/7zn1tjenp65IgRI+Rll13mfC4nOjo6nCQmpZS/+c1v5G233Sb7+vqscUuWLJHjx4+3EtmRI0cqA8iVHicefvhhCUAuX77cGtfX1+ecUKgolUpy1qxZ8vzzz2dJLwlwEHFSOHHihCyVStaYlStXSgDy+9//vjXuO9/5jmxubpbbtm2zxv3VX/2VnDdvnnMy/tOf/uRMa6jARu5Cli/Y1BKLFi2S6q2DQwGqTrhe06ku9rh80qGMjo4OZDIZZxm+853voLm5GX/9139dpZz1t9dLL72EefPmVfVLEah1kiI62tvbMW7cOJZ2lTX6RqQkIYR4Xkq5KPSzlNxTpEiRYnjCRu51cVomRYoUKVIMREruKVKkSFGHSMk9RYoUKeoQKbmnSJEiRR0iJfcUKVKkqEOk5J4iRYoUdYiU3FOkSJGiDpGSe4oUKVLUIYbEJSYhRDuA12IkMR7AAabscCLNlz+Gat7SfPljqOatnvJ1qpQy9A2HQ4Lc40IIsdZ0S6uWSPPlj6GatzRf/hiqeXuz5Cu1ZVKkSJGiDpGSe4oUKVLUIeqF3L9b6wwYkObLH0M1b2m+/DFU8/amyFddeO4pUqRIkWIg6kW5p0iRIkUKDSm5p0iRIkUdYtiQuxDiSiHEFiHEdiHEspDPhRDiG+XP1wshFlQpX9OFEL8TQmwSQmwUQtwYErNUCHFECPHH8p/bq5S3HUKIDeVnDvo2lFrUmRBirlYPfxRCHBVCfDYQU7X6EkLcL4TYL4R4SfvdWCHE40KIbeW/xxj+r7VPJpCvfxZCbC631c+FEKcY/q+13RPK2z8IIfZobfYuw/+tdp39WMvTDiHEHw3/N7E6M3FE4v3M9P17Q+kPgCyAlwHMBlAAsA7A2YGYdwF4FIAAsBjAM1XK22QAC8o/jwKwNSRvSwH8qgb1tgPAeMvnNamzQLu+jv6LGDWpLwCXAFgA4CXtd/8EYFn552UAvmrIu7VPJpCv/wdArvzzV8PyRWn3hPL2DwC+SGjvqtZZ4POvAbi92nVm4oik+9lwUe4XAtgupXxFStkD4EcA3hOIeQ8A9U26awCcIoSYnHTGpJRtUsoXyj8fA7AJwNSkn8uEmtSZhssBvCyljHM7ORaklCsBHAz8+j0AHiz//CCAa0L+K6VPsuZLSvmYlLJY/ucaANO4nucDQ51RUPU6UxD9X576lwB+yPU8KiwckWg/Gy7kPhXALu3fuzGYQCkxiUIIMRPA+QCeCfl4iRBinRDiUSHEvCplSQJ4TAjxvBDiEyGf17rOroV5sNWivhQmSSnbgP6BCWBiSEyt6+5v0L/qCoOr3ZPC/y5bRvcbLIZa1tnbAeyTUm4zfF6VOgtwRKL9bLiQe9hXlgfPcFJiEoMQYiSAnwH4rJTyaODjF9BvPZwH4JsAflGlbL1VSrkAwFUAbhBCXBL4vGZ1JoQoAPgzAP835ONa1ZcPall3twAoAnjYEOJq9yTwbQBzALwFQBv6LZAgajlGPwC7ak+8zhwcYfxvIb8j1dlwIffdAKZr/54GYG+EmEQghMijv9EellL+R/BzKeVRKWVH+edfA8gLIcYnnS8p5d7y3/sB/Bz9SzwdNasz9A+iF6SU+4If1Kq+NOxT9lT57/0hMTWpOyHEdQD+J4APyrIpGwSh3dkhpdwnpeyTUpYA/JvhmbWqsxyAPwfwY1NM0nVm4IhE+9lwIffnAJwuhJhVVnzXAngkEPMIgA+XT4AsBnBELXmSRNnL+x6ATVLK/2OIaS3HQQhxIfrr/Y2E89UshBilfkb/ZtxLgbCa1FkZRiVVi/oK4BEA15V/vg7AL0NiKH2SFUKIKwHcBODPpJQnDDGUdk8ib/pezXsNz6x6nZXxPwBsllLuDvsw6TqzcESy/SyJ3eGEdpzfhf5d5pcB3FL+3acAfKr8swDwrfLnGwAsqlK+3ob+ZdJ6AH8s/3lXIG//G8BG9O90rwFwcRXyNbv8vHXlZw+lOhuBfrIerf2uJvWF/gmmDUAv+lXSxwCMA7ACwLby32PLsVMA/NrWJxPO13b0+6+qn30nmC9Tu1chbw+V+9B69JPP5KFQZ+XfP6D6lhZbtTqzcESi/Sx9/UCKFClS1CGGiy2TIkWKFCk8kJJ7ihQpUtQhUnJPkSJFijpESu4pUqRIUYdIyT1FihQp6hApuadIkSJFHSIl9xQpUqSoQ/z/M/qUee0SMGwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "t = np.arange(0,(len(abp)/fs),1.0/fs)\n",
    "line1 = plt.plot(t, abp, color = 'black', label='ABP')\n",
    "line2 = plt.plot(t[0]+((pks-1)/fs), abp[pks], \".\", color = 'red', label='peaks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Detect pulse onsets and peaks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Import the functions required to detect fiducial points by running the cell containing the required functions below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Detect fiducial points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "fidp = fiducial_points(abp,pks,temp_fs,vis = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plot the ABP signal and the pulse peaks and onsets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7faecc631760>]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABfAElEQVR4nO29eZRdR30tvOtO3a2h2xpaas2TbcmWLdsajBSC7aDwLJLnYPIleWaFh8Ng4MWPOPkIQV4G8gDbCEJCAgsTOwsHowcYPhKCQ2JiI7BlrMEWtiVZSLLmsS21rFnq+db3x711VX36VNWvzvmde7uvz15LS62+P9WpcdeuXVXnCiklUqRIkSJFfSFT6wykSJEiRQp+pOSeIkWKFHWIlNxTpEiRog6RknuKFClS1CFSck+RIkWKOkSu1hkAgPHjx8uZM2fWOhspUqRIMazwq1/96oSUsjXssyFB7jNnzsSmTZtqnY0UKVKkGFYQQhwwfZbaMilSpEhRh0jJPUWKFCnqECm5p0iRIkUdIiX3FClSpKhDpOSeIkWKFHWIlNxTpEiRog6RknuKFClS1CFSch+mWLduHbZu3VrrbKRgxtmzZ3H8+PFaZyNFHaCuyH3t2rU4c+ZMrbMRC6tWrcK73/1uZ9xb3/pWLFiwoAo5SsGBYrGICxcuOOPmzp2LiRMnViFHyeH++++HEKLW2XjTQwyFL+tYvHixjHtDtb+/H7lcDnPmzMHu3buZclZ9qEHhahdqXIqhgXvuuQdf/epX0dPTg3w+b4yrh3ZVZejr60M2m61xbuobQohfSSkXh31WN8q9WCwCAPbs2VPjnKRIMRiPPfYYAODixYs1zknyyOVKbzWhrFRSJIe6IXcfpfPFL34R69atSzA3g7Ft2zbcd999w1qRpYgOpda7u7trnJPkMWLECABvDnI/deoUtmzZUutshOJNSe4rV67EW9/6Vmfc97//fZw7d84aUywWsXbtWmdaK1aswIMPPogTJ06Q85mihOPHj2Pbtm21zkYsFAoFAEBnZ2eNc5I8qOT+vve9D7/xG79RjSwlhqlTp+K6666rdTZC8aYkdwq2bNmCO+64Ax/5yEescd/4xjdw880348c//rE1TuWPqtz6+/tpGXUgk8ngwQcfZEmLG6dPn8aRI0eccbNnz8Y111xThRwlB6XcqeTe29tr/bynp8cZww0pJclWopL76tWrsX79epa81QpUm+3HP/4xDh8+nHBuBiIldwPUIHR5+OrzXbt2WeMaGhoA0MndFbcUwEoAcAwOKSXuu+8+0jOrjY9+9KOYOnWqk/Coy/sPfehD+Iu/+AuOrJFx5swZ3Hvvvejp6bHGUcldtWvPs89a4xoaGkinpc6fP+8UClJKPPzww06iWr16NUaOHImdO3da45ahVAaxcaMzf0MVnZ2dWLNmDTneNdHefvvt1Vf4Usqa/1m0aJGMiwsXLkgAslQkM/r6+khxmzZtkgDkwoULrXH33nuvBCDvv/9+a9xVV10lAcitW7da45YCciUgz/z0p+agdevkBUD2ArLY1CTlunXGUEpZi8Wi/NrXvibPnTtnjeOGytvhw4dJcdT0qol77rlHApDf/OY3rXFXX321XArIXR/4gLm99HZtbGRpVwDy/e9/vzXu6aeflgDkXXfdZY17z3veIwHI1atXm4PWrZOdmYzsBWRfQ0PsMtQKH/zgByUAuW3bNmucKsOZM2dIcdwAsEkaeLVulDsV1KWs2vF3xVMVeWNjI5YCuOwb3zCr7fXrsQbA5wGMuv12c9wzz6CA0jetFDs7gWeesT7bhWeeeQYf+9jHcPfdd8dKJwqWAij87d86VyBAdY8HSinx1a9+FadOnbLGKVXs2ptZ3NuLNQBmPvoosHx5eHm1dpU9PdZ2pazcVH398z//szVv+XweSwEseuopa3pjxozBUgBzfvADa9/MF4vIARC9vaQy9D33nDV/H//4x51l6O3txbhx4/Dd737XGtfX14cvf/nL6Orqssbt2LEDAEh7ZEsBiFWrjHVSzX6ro27InVqBVHLXz+raQCX3RT09WANgyj/+I21wd3ebB8Ytt6AHQC+AnvK/wyClJJGAOou8d+9eaxk2btwIIQQOHTpkjQNo9bwMwBoAY7/yFXOdlLEUwP6PfpQ0CbjwwAMPOL3eDRs24J577sFdd91ljWtsbATgbv//VihU2hYm4tbatatYNLarLgKKb3+7lVQo7T9p/36sAfDBAwes7XDt+fNYA+AtP/mJOe6WW9CbyaAXQJ8QpDJkb73Vmr+/+7u/wwc+8AHj50DpZu/JkyedAmX16tX4xCc+gfvvv98apzbAXf14KUp9eOSqVcY6Sck9JrjJXSkyCrkvBfC2X/7S2kFv7OxEAYAoFkmDu1tK88BYtgzLAXwGwP9z2WXAsmWhYXLdusoAsg3alpYWAKUNThtWr14NAPiXf/kXa9znP/95FAoFpzq6BUABQBYw1wkuDaBpjzzinAQo+NSnPuU8paHuTbS3t1vjFLm7vPTDl19eaVsUCuFtu2wZfm/ECHwGwPLyv0Ohr9y6uoz1Rm3/UZs2uSceAFcfP44CgIyU5rhly/C5m27CZwDc3N9PKoPtmVQogeJqB7U30tHRYY1TeySuTdBbAGedpOQeE3oFvvTSS8Y4H3JfCuBDHR1WMpnZ3o41AN65fr11AO2dNg09AIqZjHVwv7etzT24ARycPBmrADSYJgCAPIDUKsVF7rNnz3YvywH84Ac/AAC8/PLL1vTWZjJuwsOlAeQqB3WTmQKl3Kh2W+aLX7Q+d9f48VgO4LPZLOTPfmZs24ZbbsEqAP1LlpgfqomAXos6Fs8+S6q3c4sWVdKTlnY4uWABegD0C2Ftr0NTp2IVgA3mEtCFDBFKjLnai9quyqra8Sd/Ym3X3VOmoAdAXynx0HLo3FTNM/FOchdCPCqEOC6EeDXw+48JIXYKIbYJIb6k/f5eIcTu8me3JpHpMOgV+JOf/MQYp5P72bNnjXFNr7yCNQD+39OnraQ9fd++kvq0qRkAJ664AssBPNTWBqxZYxzcHZdfjlUAhIXYgUt2kG2yKt50E4k8Vd25yH3uyZOkiWzChAlYCmDsww9bB0bzrbdiOYC/HzPGWifPAO5yaMt8W96oKkopN9cpmNnHjmENgHs7O63P7evrwwYAD/T3o89G3GVYVz3ayu2+t7zFWG96+xfzeWP7n7/22kp67atXG9M7fdVVWA7gU1Ja22v06NHmvIeU4beKRauQoUzaxWKRFKdWWurGsAnzTp0i9aeto0ZhOYDH5swx1olaBQKluzPVAkW5fwvACv0XQojfAvAuAAuklPMBfLn8+6sB3AFgfvn/PCSEqMrLJfRBazs6p9sststHI4lL1aNXXEEi0BEjRmADgI8dPWrtyEqBuIhWldemQOTSpZUB9No3vmG2b8ppnT9/3vrMya+9RtrwW9jdjTUALv/2t502ygYAX29uttbJgUmTsBzAozNnmkmFuBmpHwu0WW5q5fa7W7dCWm4zX374sPdzrW1WbguXpdUxZw5WoSQGTCi+5S2V9v/+XXdZ238DgFUAjkyfbs2biuu/8UZjnP7SMFsdbxTCrfC1SVta+lJm40YSGSsbdSWAzp//3PjYm6UkjX9VJ49OnOgcXwCq+kI1J7lLKdcCOBn49f8CsEpK2V2OUe8ofReAx6WU3VLKfQB2AzD3AkZQyV1Xura4U+UlqIu0X581qzKAen/6U7OKKs/es2fPNj4TuEQCrsFNIfdisVgZjP9u2fXXN97O/PSnxrjX5827VCcWJbjo3DnSaoZKZLlcDhsA/NO4ceZJQFvm92Uyxrz19fVVyvrkZz5jfGbjyy+TSOUIcXLXSY5C7i7bQMXZPGa9/Y9Z+p0+dmxetK5Ajx49SkrP1ra3lOtq4cKFxhjq/kL++edJk+yEPXsq7dr4u79rbNd9M2Z4rXpd7aDAdTmRgqie+5UA3iaE2CiEeFYIodaZUwDoRykOl383CEKIDwshNgkhNrk2NyiIQu62Sxun5s2rkPbr3/mOkVT6+/srA+js/PnG9FQDuy6KUL1DlZ7NNtDrxPYmwhGbN1c6/Oh3v9vY4Y/PmVOpkwOPPmqsk73Tp5M8XMrA0OOsdbdsGX6nUMBnAKy+807zJPv885WyvvPLXzaWVV+52Y70HZk+/VI/+b//19pPFDjJndr+NgtSj7PVMXWM6XGUslqtL+L+Qudb3nLJgspmjXFjt2y5pMgt7bq7tbXSrsWnn3aqctskRq0PbkQl9xyAMSjZYZ8A8ANRWm+ErTlCTU4p5SNSysVSysWtra0RszEgvcrPHOSuk7ZN9eiD1nbWmXoeWpE2dXBTBg8AnDwZXHxdgn5aQhCWoKtQ8l9N0CeB3iefjDUw9DjXxLgxk8EqAK+NG2eM0TcZM/39xrLqK7diLmckC72f2CwSvZ/YyptU+3PHUcmMEmctq+bNf/M97zH2pc7rr6/E/fKznzXGHbvqqkvCw7L67OnpqbRr5/XXO8tAVe6uvs6JqOR+GMC/li9JvQCgCGB8+ffTtLipAMzrN0ZQyV0fZC5yV7B9AYjecDbiVnEXLlywLs2oyp1zcJ9duJBEZHpZbekVCoXKwOi64QZn/np6eqx1ouJcJ53U57Z27V62rFLWfovCOzl3boUsfvXFL8ZW5EnZMlTlzk3anJOFa9N674QJWAVgx5gxxhjdgto/aZIxruPyyyvtuu0f/sHYrnqeKKsUG7kPN+X+bwDeDgBCiCtROql2AsATAO4QQjQIIWYBuALACwz59IKrsyjYGkQfjDZypw5u3zhucrfViX5a4plPf9q5b+DzXE6ycBGZqjvbYNQV3nc/8AFjWdXpFtcmo95PbPmrtS0zbJU7LvU7ysoYsI9XvV1tK3JdSFDqmKrcq0nuOVeAEOJ7KB01Hi+EOAzgrwE8CuDR8vHIHgB3lt9zsE0I8QMAv0bp6OfdUsqq7CBE6aAcyl2P8xnc6s15QaiO0Nvbi2KxiEwmfP7lHtwbUDq1cPWMGSzp+cZ1dnYa64RSVqrdpgb3BgD3jh1rjVOgrgQ5lXuxWERfX1/lNRimuFqRO2ccx0RGtT6o49VXyHBMdtxwkruU8j2Gj95riH8AwANxMhUFVJVKjaMOxiidhToJ9PT0VM7lBsGp3KlxUZQ7pxKk5o16FJZKxhyKPGrduci97j13YlyUSZZjTAStxbCvFayVcq/LG6ocjRZFkXOQAOfynbtOak0CNs89yqTNMblHIYtqWlWUtIay3UaNo5aBOl6p/UmHqU5qpdzf1ORuIwtuco+q3E1IalleK+VO2ZDq7+83brwmSe5cXjrl6rsvIXNP2twbqpT0+vr6BvQtUxwHaevtyiEWKH04Ve4xkaQtw03anCqqu7t7QJnCYlx5oz4zSYVHjTMNSCrxcCt3H3JvamoixVHTA4a2jcb9XI6yUttfR9yJcbidlhlySNKCqKZyLxaLFd+O0uGllMYr3rWqk6SUO0Aj92raLT4bpZTXA3NuRtZ6oxSgnemnpsd5Gg3g68Ou78dNbZmYiNLxOJZl3J57FIVnem4U5catBCmvFnDFUfLHvSJLgiySIPdqnpahnkjh7ndJ2TJcQka1q6lOUlsmJlRj5PN59uVWtcmdQgKUDpOk0uaKcw0MBbWaoZA7B6Ho4JjwAJDK6kuMQ7VdlZqt1sZrrTZU9T5M6Zupco8AVYGNjY0sG6o6uOwWSpyPwlPvi6GQ+1BU7pSBoeJcrziupbXkagcVNxyUO9ceCef+AuWVDFH2tLj6ieqbrhW0ECJV7nHQ0NBQVX9Z7yzWd6t7ECOVBHwUQ628dFeca2BQ44ZDGXzb1eeyky2GmpZPHHXSjmu36XFc4okaRz3d5PqqTfXMpqamlNyjQFVgoVBAT08P6QQJ91GouASl4qhq1kUC3ISXhHJPgtyrffY/m80il8vFJgEV56vwqz3hZbNZclmpY6xaG6q1GK9q3DQ1NaW2TBSoRlMN4lIzmUymZqTtmr2pJEBVDJS0FCjqqFAoONNTG6Vcqpf6hcWUtIASQVHiGhoanGUQQpDqJJvNkp6bBLlT0qJMUEDJ+nRN2qq9XGWlWlqutHytD0o7cAuPxsbGVLlHQZDcXRWtFD4lPdfGmyIy7s4SVx0lZcu4OqnPsty3rK52pdpy1DjKHo4QwmsS4CJ3V79TZcjn86QVmcs2oNoLUkpkMhnkcjnWiYyi3Juamtja33c/yJQ/Vb+NjY0oFotV+8KOuiV3l5ptaGggLRkpnYB70FLjXGqW23P1IQHqspxbHSmSddlyFDKmxilyp/YTV51QT0tR+zplMlZxFI+cUgbqmKC2P0And0pZKWXgsmX0+gWqdxyybsm9Wso9SkfmInfustrS0uModZLP5yGEYFNHrmW+PoCklE51xK3cXbYMgMREgMvrpZI7hfCA2ggewP76Cd/VB2VFzrWq1NvBFseNuiN36ixaS+XO1eF9CI/zZBCl7riVG9WCSsK+oZA2dSXAKQJ8ytrb20tazVDGhI9y51ilKJsH4OvrVIFie6aK83ELbHHcqDtyp5AFUH3lrueNW7mbBhBVuemkTY3jInfqJqMPkQG0gcal3KnpcatZX1JxHTLg3K+glLVYLJLPw/tYHxzKXU0olEuRPhuqtjhuvOnInao+FDgHLfUEQRI+NEW5UYmM03PlqrsoKzdOz51yWkYI4SQLnfC4vF6u1UytVm6USSAJcucqQ2rLxERSnjsnQSklwHEUUt9Qo1oVlDjK4M7n8+zkznUyyEdFcapUztMyvpdnuMidezXrmsioqxSfMvhM2lztSr1vkm6oRkSSniu1I1MIz2czkmtwUzsfpSOrCYpSVsrg5lJH3LaMgs9pGSq5u+ouk8l4TXicyp3Lc1dliEvuvn2YuinMKdpS5V4lUAd3LU4G+CzfKWRBtSq41Kz6TtcklDtlj4RaBkpcQ0OD9cyxj8LntqCo7e/T14HqCZ5isUieyLg3I2thGaoz/dQJKiV3T9RKuQP0I24+Cs9ncFOtirhqVg1aH3LnWs1Qb6gmoWa5lTuVVKgnTahljWvLJbHn4ntGnLqXZtpfAvzGK0ecvlcBpLaMN6IMWinN56GjdPhqKzyuXXpfWyaXy1VdHfn4xgDfMp/rhirgRyqUOuYmxoaGBvT39xu/8i6J1Sy1vXxWKZRx7TNBxb3Elir3mIiyeUSJ47ZlqAovl8sZj66pOOoNVU5bRil3V96Gw0RGTY/rhipAFwGu9tc31KmKkRrn6k/c1odtTERtV1sZfPMW125NPfeYiKLcAZ6lKtU39uksLutDLwPnhir1xE+1T8tw7y9w2XcAvBQedZOZUse+fjVXnfiUIW5ZffswRfBU25ahloEbbzpyV6CcOVfpcXUC9VzK8t2mjpMkMtfxS1/PnVu5c+4vAHTlbvJwfb15LuUexXPniPPZS+HYX4jiuQNucqdcTuJafaTKPSaieHMAbQnKqVKoy3eb56orKEoZfPzlvr4+q+eahHKnrhhsN1l9fU3O251UMvY5CkmxvpKY3F1xvkdc4x6ZDBIjx6TtK8ao498lxlJyj4ioMzx1WW5TbkkchaR0FtdrVfX3rwPu6+eq7myvM0jinDvlYhd16ctNeBRS4dyHUKsj14aqTxl896FsYkGRMcfkrq8EqWWohS1DGa/ZbNbY/tQJiht1R+5JqRmXcqN2eCoJUJS7S/VQOzz1qFZStkwul3OebqAu32vhQ/sqdy4lWO2yUshYLwOnLeOz2euKS+oABOXkHZAqd28koVKo6XF7qa7BrX/rDGVgKPuGShYu5cZ9zt1WJ0Fyr/ZpGQqBJkHulMmduvrk3IzW29W1muWYBHyJkVLWJMYrpQ+nyj0m1LfTcG9G2hrEZ/lGPSNOGdzUJS2V3F0efq2Vu8/ynaNdgUt1Ysuf3l4mwlNl4FTu1JWb72Y0xXOX0n6W3LesXBMUt3KnkjvFlkmVe0SoxlXkTrnEAPCdvqDaLRxxvlaFz+Yx4LaguMndNjCoZVXw9WZdz6WIBVUnAJwXgFz7C0nYcpyrWdUOrvR8+0m1bBngUl/q6+tzrj5swkNPb1jaMkKIR4UQx4UQr4Z89pdCCCmEGK/97l4hxG4hxE4hxK3cGTYhqgVB7fBxT0sA8I6rtnL3OUZWK8/dVVbq5O67mqG0PyWOY3IP2nJcF4AocfpE5uon2WzWSow+tgy3BUXt6xTlDoA0rl19kxsU5f4tACuCvxRCTAPwDgAHtd9dDeAOAPPL/+chIUSWJacO6GpGCBHbX1bwWZZTOwHlhWU+yp3blrGVw3dZHvekQfC51Sqrb51QyT2bLQ0Hm8L3mdx93mkSN85HuaujqxRy97HbOBS+z+qDw5ahchM3nOQupVwL4GTIR18B8FcA9HXNuwA8LqXsllLuA7AbwI0cGXVB7/CUpRR1hqcOWo6N0mB6lMFtG0BRbRmO8/Xcm1Eqf1zkzq3wVXoUNet6ro8tR3mnSRLn3AH33owPuXNvHnORu0t4+PRhSp1wIpLnLoT4PQBHpJSbAx9NAXBI+/fh8u/C0viwEGKTEGJTR0dHlGwMQLACqRYEB1noEwrFw+P03CnkzrWh6mtBuCbZKLaMiwSoG+q1smWo5E6d3Clx3GfEfWyZYrHoHBO2PqxWFblcDtlslmUSSEK5U/swxb7hgje5CyFGALgPwGfCPg75XWjLSikfkVIullIubm1t9c1GWHoqf6SK9lWpXIOb+7QMJ7n7bKgCtMHNtaSlLt+z2SyEEFX33KlxFHJXl9O4J3dXHGU14ztBAXYLymXf+LY/ly2jnumzoVoPyn0OgFkANgsh9gOYCuAlIUQbSkp9mhY7FcDRuJmkgDo7+i5VuQc3t3KnlJXTlvEh91wuByklyV/m2lCl1km1lbvKm4sYKW/eTKqsVHL3XaXYNoa5V26ctgxg3yilugUKQ165Sym3SiknSClnSilnokToC6WUrwN4AsAdQogGIcQsAFcAeIE1x+Z8AfC3Zaqt3F2euyoDh3L3ff0AZfNQj+PylzmOQvoOtFrYMoDfHg7Xyo3TqlJKG+BT+JwrN1+BEseW8XULAAwt5S6E+B6A9QDmCiEOCyE+aIqVUm4D8AMAvwbwUwB3SymrUhLfiuZS7vozbXFUW0alR73EkoQtw6ncAZoS5Fbu1SQ8vazV3lClWhoUknK1l75RSimDTxzVlnGNHW7PnTrx+EwC1VLuOVeAlPI9js9nBv79AIAH4mXLH0lZFVyKzCdvKj3Xu2C49heCCp9zQw2IZ1X5Krda2TKcnju3cucid99Jm5vcuU6G+W6oqoMSQohBMXre6sFzH/KgVLSvcufwIVV6VCKL6/VFVe5cG6o+JBD3TL+CL5FRL55U+7QM58qNOln4TlBcZa0muau0fC4xAeGbwtTVov7cIe25D1Vw2zIK3IPbtsnoOxh9T0vEXc2oZyah3LkGN0Cb3H1PkHBfYoqj3PWVG2XPhXr6JqnTMlwrN847HVQbzVYG35URMMQ89+ECX18rqQ1VamcJi6P6i74dPpPJIJPJsL09kku5K2SzWdKE56NSOSa8KBulXKuZoei519KWoa5SqjlBRfXcU3L3hK9VkcQlJkqcLT3uzuJr89TScwfCB1pSRJbE/QVKnE+dVMtzV5Mq1YKspi1DXaUkRe62skb13FNbxhO+hJfNZpHJZGpyzt0UF6bc9SVdWBznbj6VtJM4CgnQBtBQJXfqqSpfz93V/hw3I6N67kNxQ7Xa1pLKm4/nnip3TySlepPw3AG3LUPdyOHcUEtqQ5XzAhCVBIaicgfgdUPVFBdlcqes3Gply3DuuSRF7rZ2UM+leu6pcvdElA6fyWRYN4UocRwXhXwnKNcXLlNtGf2ZtrgkBxCX507pJ1H2UuKqXv3Wpim9WnvulJun+gRFuRTFtb/AtfpQz/Rph9RzTwhRB3fcjqw/M25cmHKP06mC5+Grffabw4dOyoLyJfdaXGIypVdrq4rLlqG8PZJaVgVOywiwr9y5BQU36pLcuZSAz3Kb41RFcBlNiUvClqm25869oco1CXCuyDgtjVoqd5dlRC1DMI7LlhFCQAj7i+OS2A9S9Ru2R6KQKveIiKrcqMsyLvvG57SM6bm63eJzdbvap2WSUO5xiUyBoraikHY1LgBFaX8Ocg++fmAoXmLyXZFzHYVU6dn2yIDUc4+FpDZoqrHxFtwPAMI7n+/rB6iDO4lX/prKoMdx+MsK1DoB3JtgXHabnh6HbeB7PJDqayd1jLCaG6pU0ca9quSwb7hRN+TuOzuqBqa8jhao/iUmn2U51YKI+24ZlVZSp2U4LAiA/ppWLltGPRNw1x0HMVLrxHcSUKKCa4LymfD6+8O/6CbKKoUSl8QlJlN6qeceE74qFUBNTsv4LvO4PFefo3Acx/mS8DXVZGwjAS7PXaEWRyGpy/xaWRpJrVLCyhplgqqFLZMq9wRBrcCkOjLXZSc9b0B1N9RcCj8p5eajeoB4L3FKQrn7tD+17my2XNT9JdsqFbj0Dppq2zLcfZ2T3H2EhylOIVXuERFFuXGSO+clJsCu8JIgbRXneqdNkoObazXDacsokuXaKOVQvVEIz7ZKDarjah2FTKqsnBMUpy0D2O9qcKMuyb0Wyp37nDv3spzrAhDXZhQ1LoqvyXnOnToJcBKe69uOhootE7es6pnVPgrL1Yej2DKpco+AqIM7jqpU4PJcg0qLEseh3Khx+jNNeVPpcalZ6tKX23NXeQNop2q4yD143HAoeu6+7Uq127hOS9XyKGTquSeAKIM7k8kYfUg9BuC5xOS7BKXEURR5rU8a1MKW4VqlANUld1/PnZvck/Dcbe3KvUoB7O9Mj7qqiGPL6Omlyj0CkrBlqOkl5aVS40zPTPqkgesYKYfCU6iV5w7UhtyTsGVcG6pUS4PDlgn2Te6ycrQDANIEpdeJay8tVe4RkJQtA/hddqFuRlKW2wCtw1OPB3INDJuq1OM4lTtAt7S4bRmX2vJRgpzk7rMio9pyrvR0pW2rO19FTo2r9mkZ6l6aGhOUcZ0q94igVmBS5E49fcGhUvUBROlUVOVWK889Dgno+eNW7lxWFUD/2j4f1ctpy9mUZXA/IE5cmOql2DdxyV2PsZXBp28CII/DVLlHAHXp4+vNqfSqsSyndpbg6wcA2usM3gyeO0Bb+nKSexSrijqRxSGLqLYcly3DMWknYS0C/pfJ4rwuRCFV7hERZVnuUrO1Ivek4lwkoNRRHIIKi6vFJjPlhjKFyFQcl1VFLYPPhiqHBUmNi9Ku1b6NzWnLcPZNIPXcIyFKJ6C8fgCgkQCH5xp1cFPjuOpE/eFW7hSVSi1rXMJTMUBtyJ17z8V2Mmw4EKNPH6YcgKD2Yd89stRzTwDc/nKQ3DnUrErLFBdVpVLT464TLoXHPZFxbagrJFEn3Cu3OOfhfd4yqtrLhxhrsaFajRVZFDGWKvcI4PaXVQxgf8GYQq1tmTgXQKKqWeobNbnrhOJ/ck9kXGRBmYx9N8pd6SVVVi67zeeiEJcto8rA3Tdd92ZS5R4BSXQCBU4f2ufIlCm9KHGcG6oAr3KvNpFR43xIwBUXLENYnegKutq2HGdZfSYyrjIocJeVOl4p7eV6Jjfqltw5lRvFh+bwXBW4CY/zhqpPnVCVexJExtn+1JuslDqu5YY69fBANWy5sFUKx6RNzZsrTn8mNW9APHuMG3VJ7mpZ7tpkonZ4SmehqGNqZ1HPdMW5CDSpG6oqf0NVuXOtUnxIwBWXtC0Xt058bZm4J8hqfUNVlZU6GXMo95TcIyKsw4eRu4Jvh+eaBGrpubvULACnX81VVoVaWBDUOB8SoMRxtX/w6KorPcrKLQlbhkOgJLn65NxQ9Z0EbNzEBSe5CyEeFUIcF0K8qv3ub4QQO4QQW4QQPxJCXKZ9dq8QYrcQYqcQ4taE8j0I3P6XrwVBTc9HaVHiqum5qrRccZybUdSy6vnjtmV8+gnXJSaq18u1cqOWVeXLp/2p7erah+AgdxXjW4Y441V/LqUPc4Gi3L8FYEXgd08DuEZKuQDAawDuBQAhxNUA7gAwv/x/HhJCZNlySwCX/8Wt3IN54yBtH7+Sa4JS8NlkjrOkVfCxtDKZDKSUJFtuqKrZJFZucfcX1GsFXHFJ7qVwKHcFHxuV45w7YLd5uOEkdynlWgAnA797SkqpetMGAFPLP78LwONSym4p5T4AuwHcyJhfWz4rP3O9LpVbuasYat5McUn4lUkQGefgBlDV0xfDjdy5NlSpnjslLqkyuCZtztdKcJ9zpwgZLnB47h8A8GT55ykADmmfHS7/bhCEEB8WQmwSQmzq6OiInYnhYlX4HoWsxTnnWhEZ9/FQypnjateJ+jOUN1Spk4AtLqm3QgK8K/Kh1Ie5EYvchRD3AegD8B31q5Cw0J0DKeUjUsrFUsrFra2tcbKh0lN58jqb6iIAwN2RFTg3Gbm9WVsZisWil5cK+B0P5TxpAMR7iVOwvVxvZ1Rx3CJgqE3uqs6ok4DrO1m5LQ2fFV4S5E79Uh/quK6GLZOL+h+FEHcC+O8AlstLNXsYwDQtbCqAo9GzR0dSHh5Q/Q1VlRa1DNS4aq9SAHPdVUO5u+J83h7a1dVljPPdr3DVSSaTqapfrXvptbDlOMqqYmpRBpUWJW7I2zJCiBUAPgng96SUF7WPngBwhxCiQQgxC8AVAF6In003kl6+cV0A4dh9T3JZTrm6D8QfQEmWtRbLd464qHsp1Lg3gwXFXYYkVtCm9LjhVO5CiO8BuAXAeCHEYQB/jdLpmAYAT5craoOU8qNSym1CiB8A+DVKds3dUsrkS4HoZFHNDdWkBi1nWZOIq7YFwd3+rtWMAueExz2R6eerVT588hYnrlorsiTKUIv9IC44yV1K+Z6QX3/TEv8AgAfiZCoKfCuwVoRXC5Wa1CUWymrG9dK1elLucS+A+ZY1qn2jTh3FLUOc9g8rA6UP28rKOV6BaH2Tw5vnQl3eUPXxyTg7MqWzqPxxqVlbWYPLfK6bpxxxSe4vVFO5J6VmbYTna98kQYxxV7PUMviUlbu91HMp1hKHN8+NuiR3H2Kspi2j4LPJyBknpflij88E5YpLypZxvYgsmN5QInf1TKB67V9LYkzSgopTBpWWbxm4rcUhcYlpuCEp5VYt1avnjVoGSpyrTnwnqLgKL6rSdsWpvAFDZ0NVfyY1rtp7M9SyRjkyW609F/VMn69ZrHYfTm2ZCIi6zI+jtIJx1bYqOC+KVHtZnsTgVnFD0ZbRUQtyd/naFCHDeWQyicl9KKw+UuWeAJLYUFPgVrNxO3zYlzrEVUdRNtS4/WXqAOL8wmWARni1Uu7VPmnCdWR2qNoyUYUM1/uRUuUeAUkrt2rYMkN9Wa7ScsVFXdJSTjfZyhCM41rm14LckzpBVAvVO5yPQtrifPucSstUBm7UJbn7bKjGPc43VJZ5cQcG1UtViDvhBQnPNQmotGxlUOlxvX6Va8ILK2vcY38qPY7LM1FeP8E5uXP0dRXjkzefF4xxrz5TW8YDw0m5+2wyVWNDNanr5yrOVVZbeknViXomNY4y4VHa1ZZe0jdU9WcE42q5Sk1iQ32oTlCpLRMBUWbHam0eKiTZ4YeKBRGMowwMlV41bgEOF1um2pvMUfq6662QtjjfMlAuO9XSlhmuX9YxLOBLjECp0aSMd/bbJ05/rqsTqHJUaxKo5cAAeM9+D7XTF0l57lT7ptqXmPSVYFzVO1TO6nPsB6m09HIlibokd59OQInzORkQ50sCgpMMl1XBNbjVM1Wa3LZMnEtMSdgynGSh8maL87HRVJzLS6fE1XJy53wDZnBcxxEy6pkqj66+xPF6a27UJbn7LN9McSotFTcc1WyUwc39VkiqLRO3rAq+qxRTekm0q0Jcbz6JfYhq9/Uk9hco7arSAuJfduJuB27UDbkr+BKZLU6BsixTcbXyZrk811pNZFzeLOBu11qRey3bnzrhcd9zqObhAWq7KsS1Zaj1q0CN40LdkDvV14qyfHNZEPpzuZblKi7OC6aSVm5xXrqWpOfObcvVitw5idFnwisWzd9T6tooDUuvFhvg1Diu1SdVKKa2TARE6ci+Hd703CQGLTU9buXGrdx9TstQy0C9oUqd3E3p1Zrcfbx5bmI0bW5SBE/UMgznSZtjguJGXZI7t+fK+R4VjjhVrmq/W4aSN2p6UVYp6g+nSqXGUTfeuCZtatxQW80kZbdRiFH30uOUISyO6xIT5TY2F+qS3H0HdzXUrMobR1zSg1tK8/FQBc6BASSj8DmJDKBdLecgPNdt0STaX3+mraxRVqnVWn0oUNuVekOVUgb1J86RSW7UJbm7Ni1qZcsoUOOqdTOu1n6lSo9L9VJPS/koPGpcHMLTT5Co9Lg9d4rdQi0r5yrF1U98j/NS4zhXnz6v0EiVuweibG5wbagOBc+VS7lRiZG7TqiTgOlLyKMs3zniOFduSRCeSo+bGKs1aUc9q0+N41x9xrVvuFGX5O7TQSlx3BuqtfQrXcrNZ2BwnpZxEaMC1+aWSgtI/rKT78qNk/CA6qre4Wa3ca4+ue5qcKGuyT3OhqpKC0hGuVfztAzX4KbkLZiez8DgtmW4PXfOTcZqEF4c1cvluXOOCd/LidVeafmIsdSW8UCcwR11o0xhKKgZzvdy2OIUuP3KahHeUCF37smd6yhk3MtuvqsP3Q7kblfTZqlu8Zli9Diq8KAcD01tmQiIoty5bBmFoey51/LWJpdf6XquAsfKbShYFRRS4TiCF7Ws1BM/1VqlcgkUbhuNavFxo27IXSEJIqvWhqpeBmp6HJ77cPEr1XNdl5iSsmVcl53eTOfcqSd+qvnK51qvtHxuFKe2jAfCZkeu1w9ks/FfDQz4zfDAm+9K/lBTeFz7EPVE7j5n8H3ikrBlkl5pKXCttLlRd+Tuo9x84+JuMilwEp76wzUwXDfofFczcf3KJHxobltGIYkNOg5yp5SB0v7BU1VDadIOpkUpqxJtwfxFsVtsZdXzl9oyERCnEyRtVSSl3Hx8zSSIrFbHQ2uxv0CNi0MC1WjXavd1rkmbeqbfpwymS2xxbBlq3lJbxgM+m0cK1XpNQRRvzpZeFF9zqNsy1DqJe4mpluSehF/tq8gpdgs1Pe79hWqd+HLFxZmg0ktMCYB7cKu0gOodD0sqTi8D92kJ1z5ELTcZa3FaRsrBezNJrtw4VilRL7HVg+ceFlerMnCjrsk97oaqArdyr2YcJ5EFnwnw1slQusSkPzPuO2i421VXs+qPSfX6KnKfdqVMZNV8vfVQbS8Fl6vADSe5CyEeFUIcF0K8qv1urBDiaSHErvLfY7TP7hVC7BZC7BRC3JpUxoOIMrg54tQzAbOajdIJbHE+xEgldyoJKAxVb5aijqh1olArJViLyT3uxFiNFVm1lbvpslMwLp/PW+NceeMGRbl/C8CKwO9WAlgjpbwCwJryvyGEuBrAHQDml//PQ0KILFtuLUhCufuSgFIMXLvvHCdNKKRNXZZHUYJD0ZZR4CYLEzHGWebHrTvuvu5KL4q1RC0r9+pTjdcgIYeRdm9vrzVvKj1X3JDaUJVSrgVwMvDrdwF4rPzzYwBu137/uJSyW0q5D8BuADfyZNWZTwDJKPckLI1qKLcoXqrPWyEp6VHLUCgUQgeGgp6e7cq4njfKNyxRy+AT5yK8ak54SfV1n4nMtVFKLSvXvpECdbxSyZ0SNxw2VCdKKdsBoPz3hPLvpwA4pMUdLv9uEIQQHxZCbBJCbOro6IiYjdB0jTNyOb/eHZnLw1OoxbI8iQ1VanqUMuTzefT09BjjFCjp5fN5ANUhd70MQ21yT7pdw9Krld2WhC2j4KPcqV8BOZTJ3QQR8rvBxykASCkfkVIullIubm1tjf1gvUHU4DY1SFRyj5reUFBu1RoYwThqWakDiFInhUIBAL294k4CCkN5ch8OnnscWyaYFmcZOJX7kLJlDDgmhJgEAOW/j5d/fxjANC1uKoCj0bNHR5hysy3zAbPnptLTG82UXtIdvlobqknGUZU7VR1RlbtpJVDtOkmS8FzpDYWyclwASuISUxLk7lLuNqHIjajk/gSAO8s/3wngx9rv7xBCNAghZgG4AsAL8bJIA7UCqaSt0tLjoiq8qIOb2qlsA4N7cFM2VCllDcaZPPcodce9ckuCLDo7O3Ho0KHQOJ/X4KrnVsOWC9tQD6YX5VW+tjIkeYnJZ0OVclrGtaGayWQqq8ru7u5BcdygHIX8HoD1AOYKIQ4LIT4IYBWAdwghdgF4R/nfkFJuA/ADAL8G8FMAd0spkzeXMLCiFcHHIXddCXDbMlRyHzFiBLZu3Yqf//znoXFqANkmAU4io2zQhpXVdTxMleHChQuD1HaUulNEUCty37x5s7UMI0aMAABMnz59QFxwkzGfz+PUqVOD6ipsEqjWUViXfRMsa0NDg3Uvxaddq3UAIozcL1686OybpnEYvJdQKBRC64QblNMy75FSTpJS5qWUU6WU35RSviGlXC6lvKL890kt/gEp5Rwp5Vwp5ZPJZn9APgHQVG+Q3N944w1SXNLkHoSaVJYvXz7g90ESKBaL+N73voenn37amLekNlRdqsekyIMKL5/Po6OjA3PnzrWml81msWHDBnznO9+J9Nww5UYt6y9+8QtS3IoVKwbF6HkbOXLkoHTC4pqamrB7927cfvvtA+LCTpqcOnUqND3Oo7BhdXfixAlrGRoaGnDu3Dm8+OKL1jjK3owi+FqQe19fH6ZOnTooPT0ul8vhhRdewJ/92Z+FpqdQKBSGhnIfLgg2SGdnJ7785S/j29/+9qC4IGnfddddWLt2rTHO5s3rz1Sd5etf/7o1b9lsFmfOnMHOnTutcWEDNixOLe9///d/31gGX0W2Zs0aa1xjYyMA2uDu7u7GbbfdZo1TbbF//35rnJqI3/ve91rjhBD4m7/5m9D0gnVy9913D6rrsLj7778fx48fN8apugMG1kvY8l1BJ7RgnOpvTzzxhLWsx44dww9/+EM89NBDA+J00qYKlIaGBgDu/QrV/tdccw3OnTtnzJtK78YbB56Ijrqarcb+gqlvBk/1BUlbtf/XvvY14zOBS2MiadQtuSt84QtfGBQXbDQA+OUvf+mMc9k3OgmEIdip5s2bZy3DyZPB6wXhcSpfaiCFlcG3w99///04duyYMU4N7kWLFuH8+fOD8qY6uvIYf/KTnwyoqzBrIfhZWFkPHjw4KP9hcRcvXgQAfPrTn3aWFQBWr149KM2wuMOHD5PSCzveq+J0MrQR45kzZwalERanJgEbqTQ1NQEAPvnJTw7qV2Htetttt+Hs2bPOOAB4/fXXjXnT4/S+F7bSevbZZ/Hd737XWlYfcqfukf3sZz9zPlPBNhkHJ349Tu/fFy9exEMPPTRoNcONuiN3SlwYuesNaIr74Q9/aE1PT6Ozs9OYN72hw8rgq9wV9IFkyts999yDI0eOkMqwdetWY5z+LH0SCFoG+oSjCDcsTp8g9Li4dRL8t6msYXWnoMfZ6k4RKIABxBhs/9OnT1vjVHpBcjXFKegrgmDe9D7393//98Y4vR42bNhgjNPLqk8WJuXuilP//uM//mNrWX3I/c4778SFCxcGxQZXUB//+MetE4+ehj7hBq1FvV116Cso4BI3PPLII6HxXKgrcg92doBG2sBAMg7GqU7wla98hUyMNjWjD0JdXQTj/uiP/shYVj1OIajcw043AMDnPve5QempOD0NqnILU+66IlPQ9zZsRGZTs67lbLBOwnzcsPaylVXvJ2GEq+LGjBkTGhcsw4c+9KHKZzpZRFXuCqNGjTKWISy/YXE6aQfJytT++oRr6+thVpXqd11dXYPyGZZeNpsNJWzTOAzuk5gmbZuNpq/Ybe2q91tT3nSElYMTdU/uYSo5bNDqBBVMT49rb283xumdxTZodQLVGzgY96lPfQorV64EMNADNQ3uoG0UNkEF8wkMVBY6QQQ7nz4JUMldL6uNBPSTIzZyV3sBzc3Ng/IWhrA6UaCS+4QJEyq/Dw5gPT19o9RGAldeeSX+67/+yxmnE61NWSrYbDkdtn6it6uN3E2TgE316nUXXLmZJu0wEfDwww/jW9/6ljFvVNGmxx09ar6O84d/+IeVn8PaS40JvU708pjageo2REVdkbvJ7gjGKUQh92An0ONaWlpC0wt2UF1N2ogsk8mgra3NGfdP//RPAMLJOIxkbct3ndyDdaJPAi5yD1sxhE1kKu6BBx6obLraCO/tb387/uqv/sp5LC3smcGy6nmztb9+g9o2CYwePTo0LixvanKyiYAf/ehHmDZt2qD8mcpqIzIdYRNUGGkHVw6murMJFP0sf5QVWbCfKAStJT0tfRyGbYKGTdq2CWrp0qX4t3/7NwAD2zU4QT3++OOVz4Jl1dvhrW9966BnJoG6Ifegr6Wg+7eAmbTDllRhqtc2uPVOZevIuhUTJBU9DrhEGLbB/aEPfQif/OQnrQSlp2mrE52gbMSoD0b9ucEOry+3bXENDQ348z//cwD2ugNKE1BXV5fV0jpw4AAmTJhgrRO9XW1x2Wy2skluI8ZrrrkGDz/8MAA3uau+YoubNWsWPvWpTw3KXzDuhRdegBDCmjcdtjh90raRu67+bXn7kz/5k9DnUi2NYD9RsO2R6BOUrawTJ04k5Q1ARWTZJuPFixdXVhTBdtUnp2effRZve9vbjGXmQt2Qe7Ajqx1wW+O6lLuCbgHY0lu8eHHlnLats+gDIyxOh1LSFMLr6ekZkHaUwa1bCzbCu+qqq0LjwkjWVtYwNeuKo0x406dPx5IlS5wT3pVXXuksKwDcd999aG5utk7uAPD+978fQDzlrhMBpf2XLFmCO++801mGJ598clBawbhMJlOpE9vk/pu/+ZuVk2i2dliyZAl2797trBN9s9V2IkVBJ/CwsirYytra2oolS5YMirO1l82WMcUFhWc2m8XYsWNTcqciODsuX74cn/jEJ6yDMZvNYtasWQDsg3vmzJmV39vSa2howFNPPQXA3lluueWWymeu5bYiMgq5A4OXyGEd3laGTCaDG264YVBawbgxY8ZUjn7ZyqDIzhSnt5mtrL6Ep+JcanbHjh24/vrrncQIlAauK718Po+mpiYyuduW+YBfnbjKsGLFCixdutRZhp07d6KtrW1Q++sb9JlMBitXrkQ+n2fpw/q9BdtpKYWgTRMsgzo04Orrym6JQu4+7RXM/+jRo42nobhQN+QeZsuMHj0aXV1dVjW7e/durFixwjowMpkM9u3bB8C8fFQIIx4Fld473/nOysZgnIGhd3CluIPp+aoZAHjppZecloaprMG8zZ8/P3QSoA4MWxyl7lxlEEKQiFGlR7E+whR+MG+m9jKVlVInYWUIkmAYqYTFjRo1yjq563Ec5P7Vr34VX/rSlwDQ9hdc7fDpT38a06dPd660qOLJptwpcRShwI26IfewDmpa5usVnclkSIN75syZJOVm6yw61FXmKBZE2OBWRGsj9+PHj+Omm24iEdTIkSOdZNHY2IhMJhM5b6aycip3CmlTicxEjK6BG9b+mUwGI0eOZLGgVBl6e3sHbDSbBE+w/cPiTO1PJfdgP8lms9ayZrNZTJo0CYC9Xd/3vvcBsCtyBYqNpiZZSl8SQoQqd5ctkyr3mDB1ZMBd0ZTlu0qPsixvaGggDUaANmh94my2TGtrK0nNqPTCluUu1RuWt7BJICxODTTq0pdSJ+fPnx9ArhSCMsVRJncV51J4qhzc7U+xA6iTexTlHtZeQghSnVDK+thjj+Guu+4i9WFKWZW4c5G7KoNtQxUwk3uY8Ozt7U30NQR1Q+6mQQa4Ozz3sjwY5zMYTXFUYnTZMj6Dm1InwTjTwDANIL3TZzIZ9rorFosDTuxwt3/wmSqOQu4UUvFd9XC1P6ctE/bcOO1KUeSmuKh5C0svrA9TNlTVM4Nx3Kgrcg/zF4HBFehSUCq9JMk9jIzD8ue7oRocaGF1Qh3cUSyNsA4flp7piJup7lwbrwq+VpWprGF5M5FKEJzKndr+PqtUymomqnKPQ+5qTNjOzQOl+r1w4cKgi11R2ysquYf14ZEjRw6yb6jCkxt1Q+5hs6PtaJ2OML8ScHdQlV6UjpzNZtHU1BQ6MHQ0NTVBCBHJljEphp6enkE3XqMO7mCcibSjkoCvhx98JhB95RZE1A3VOMqduvFKXaWOHj0afX19pPaP47lT2zXsUEBYf3Kp47BnUhU+dUVGaVchROj4T5V7DFAr0NRBAfckwEnu6rlUS4PTlgmLi3paIkjuccqq8hdFuVNXM6a8USY8lTeX6uVU7jYREFW5A+5JIAnlTvXcKco9rKxBxB2vQbS0tDjvJaj8UTx3IFXuJFAr0NbhOf1KKrlzEaOPBQEMPlpHVW7BOo5L7i7bKEy5q28yiupDB0G9I9Dc3Iy+vr5BL7kykXvwWVGJkRLno9wBWj/hJHfTCSKqLRNMC6AfNwxeiuLy3E2r1LA4Uzuk5E5AnNMycTajgGiKDPBTPVzk7jO4be+WUTCRexTSBmibVrYTDsEyUMtKiQsOSBsJ6CchTEqQssmsnsul3KmXbFS7uoixGp57MC316oagijbZsr7jmkruPnGmvKW2DAFxOnJc5R6EaRIIglO5FQqFQa9DjTu4L1686BzcQYVfDc/dlF7UdvUldzUgbYObEhckbVtZXX61qV3D7DZKWRXRBr+XIKwM+iQQh9zVioyyggJoyp0Sx3laxhSXKvcYCOvIuVwu9Cp4nMHt483qMUC8ZTnlHC434UkpnYM7rrXkUvg+ZJE0uQetDxe5u+JUu6rP4/STMOVOtQNs5B7c3AybLPR+4movW1kzmQyamprYlburrFRF3tLSMuCUjs/qM1XuMWB6K2SY1xd1k8nnvRFRBq1CWKeibORQyD3O4LYt34N5c5U1zlFIallracsANOWun8OPQ+5KyESZ8EykDQx+v5ArPVud9Pf3OyeBsP4UR5FT47q7uyuvKQlbGYU910d4mOo3Ve4EhDUaEL5LH3Wp6kPuxWKR1JFdm30AMG7cuAHfYmQiRsrA4LaqgvaNi4xdE54aaGp15DMJRJ20fZU75QieHmcjbT1/cfZc1HO5BErw9JWrDK446sRImbSTsGXC8kYld5ctEyY81esnUuVOQBhpA+EVbbrs5DoyR70URR20VFtm7Nix1q+os6UXh7R94tREZvONKROeydKgKPdgTFNTU+hrD7iVexBUbz7Y7+LYcio9LsILrtyqRe4UgTJy5EhkMhmyLUO5xKjnjarcbbaMfkrHJDyTfnlY3ZC7yZYJ26CjKvcgqF4qddCqd1W41Oy4ceNw8eLFSMv3uMo9ONBMEyNlIgMG14krPapyD2tX6j6E7x0BTluGEhfnuGEw5rLLLgMw+MuqTXHq24JqpdxN7erja1Pur+hxcW2Z5uZmSCkHTIwUbuJG3ZC7bXZ0KXe1S++7pI3bkceNG4eenh6nOho3bhyASwPSpBjDbJ6oalaVIfht72GrCuDS96NS68RG2nqcbRJwkQAQfTUTljefwa3HmdKjEijVlqEo91wuh5aWlkErQVOfU3FxJ22qpTFmzJhB37VrGtemLxBXiHo81ETuwW/P8ikrhZu4UVfkbrJlXApP+V9cys2H3AH3AFIEqseFdZYw1ROMU2rWNZGNHz9+wDNNcSpvwYknWMc+RAa4J4Gw44FUcg+7iAUM3mQ0Kfyoyj2IMWPGAHBPjM3Nzejp6XGu3CjkDgzew7G1/4kTJ0hlVURLnbRNcWPHjnWuKoAS0XJuqFLyRi1r2Ao/rG+myp0Imy3j2txQcVEsDcDdkU1xVHIPiwsrA1X1UDYj1TPV4DbFBSceU4c3kUVUm0e1l/rc1K6UZb56TXPw2J/+JdAqr3p/MuWtqakJ2WzWGUcld1M/CWLMmDGDvujZRO4UWyaTyUQSHmFx1HaNqtzD4gqFAhobG9mUe9SVm+0kX6rcCaBuWsRZvvtuqLnigoNWISq5K9Wje/hUcg9bHgshIit3F7n7qB4gfEM16GtS2pXa/qY4ve5MZQ16wqY4tZoJknvwucGJ1lTW8ePHOydjlZ6rXTOZDMaOHetU7lSBQiX3sWPH4tSpU87NyJaWFtJEFuUEUVxrKWwvJfXcY8C29NGP1tlIgMtzNw3aqIoszHM3Ddq+vr4B/idF4Ycpi2w2izFjxpCVO5XcqSTgOpEQrDsbaUeZ3OOQOzBQlakyBFcCo0ePRiaTqbSFa9WjlzXsmePHj8f58+cH2DdhZQg7feWaBExlLRQKGDVqFCu5SykrqtwlZFxloFhVpj4XJiiEEKTTMnp6cfYN4qBuyN229AEGLrnCOnzYcjAYVygUUCgUnIpMDUb13aHV8typcVSFN378+EF5C1Mz2WyW5LlnMhmnEqSqqNbWVgBAR0cHgPikTY2jKHIVp8qgbjSG7fVcdtllZBGg152pvfQ405igKHeVnqu9gum5Nl5ddRcmFuKUgXKgQtljQWvRdEpHn3jCyhBG7qZJ9vTp0wNe8cGJWOQuhPgLIcQ2IcSrQojvCSEahRBjhRBPCyF2lf8ew5VZG2yNCwzcBAmr6AkTJlTIGCgNSNPSl6JmxowZ4yR3ql85YsQINDY2ksndNTCiLN9tFoS+fDepmWw2G7rMj3oUkkruOnna4lpaWkhxYeTuijMpd2CgqPBZ9YS1l6oT1yQwbtw4nDlzBn19fZW4sDKEtb8pzjUJBF/2FlfIjBs3DufPnx+wIqeSe7Adcrkcxo0bVxmvJnIPpke1b2yTbLFYTEy9RyZ3IcQUAH8GYLGU8hoAWQB3AFgJYI2U8goAa8r/Thy22RG4RHimQdva2jqA3MM6ATBwErCpmYkTJzrj8vk8mpubnQQKDB5opkGml9VF7i5vnqrcqHWip2ci7Xw+j8bGRicJhKlU06R94sSJinq2Tdr6hGeK0y0NKgnY4ijk7rMiA2jkDmCAHRS1vVR6lD5M2YwOTtquMrj6OnXPbeLEiTh27JizrGHtSll92soQ3HPjQlxbJgegSQiRAzACwFEA7wLwWPnzxwDcHvMZJNg6KEAjgdOnTw+49m6KC3aCOJMAdWDoHqOprD6qp6urCxcvXnTGUfIWhdxdJMBly0ycOBHFYnFA+4e1V2tr6wByt8WpZ5rsFqBEAmrDzxZHIfdCoYDRo0d72zJUUqGIAFe7+vZhU9zEiRMBYMAYswkZVxn0dgDMk3aYGDO1q8uWCZ7Sodqo3IhM7lLKIwC+DOAggHYAZ6SUTwGYKKVsL8e0A5gQ9v+FEB8WQmwSQmxSgyUOqB3ZRtoABgxcl31jG7Q+hKcTlCm9IDHGtWUA2uDu6OhwDm5d9dgGRhRyN6mjlpYW5PN5ErkDbrLQy+qKO3XqFPr6+qyTe5gSNNkyrrP/6rlU5e5SvVQR0Nraiu7u7gEvwYsjAvSJ0RSnxqE+dqjEaHqm3ueklEYxFmyvOFaVbvNRVx/ciGPLjEFJpc8CMBnASCHEe6n/X0r5iJRysZRysVJhcWCyKnyUOzBQCZoG7fHjxyGldCp3CuG1tbVV4tRk4SILH3Kn1IkpvQkTJpAGtz6R2Za0YeRuOpHiUj1CiEETY1xy7+3tdW68t7a2QkqJkydPWkmgra0NZ8+excWLF53K3aVmgcG+tqn9hRBk5e4SAaruVH+35e306dPo6+uLTe4jRozAqFGjnH2demt7woQJOHfuHDo7O619TlfurnZ9/fXXB8RRBI/NMh5yyh3AbwPYJ6XskFL2AvhXAL8B4JgQYhIAlP8+bkmDDSarInhKwzZoAQxoYNMk0NnZOeC9zqa4kydPore31xo3ceLEQZ0ljNyDnSosrYaGBowcOZJ0ZA5wk/ukSZMAAK+//rpzYJw9exZdXV0kW0afGE2ThU7GYTF6eiouDrkHNyNdcR0dHdbJWK87V7ueOHFiQD/J5XKhZXUp7bBN6zi2jBI8x44dc07uAAasfKjkTrFIOMrQ0dFhJe0JEyZU+rCL3I8fP47+/n5rWYPjdbh57gcBLBVCjBClnC8HsB3AEwDuLMfcCeDH8bJIA7XDu5S7brm44lzLcgADNvNMg1sRhW0SaGtrw/nz53H+/Hn09/eHpgUMHkC2TuUigba2NgADCSqOwtPVsS1u+vTpOHjwYCVvYfURLKupXVUZKModcHv4+iTgIgEAaG9vt8ZNnjwZUkocO3ascnolrG2DdgBlwotLjPrEaGuvyZMnV8rqIveTJ08OsLRckzuX3aoIGTBPKCrO1a5qD8e1Im9vb7eW4bLLLoMQYujZMlLKjQB+COAlAFvLaT0CYBWAdwghdgF4R/nficNGArrqoZK77bSMinMpdxXnmgRUZ3FNAkBpoNnIXVcMcT13ndwpys1FAvpzbQNj+vTpaG9vR3d3t1W5U8i9ubkZDQ0N3srddqoGGKjcw+J05W6LU8R49OhRa/sHbRlbX3eR++jRo5HL5ciq1zVp62Wg9BO9/SmWRljMyJEjUSgUBoxr0wa4KgNFjB07dsw68VAFz6RJkyrjxiY8L7vssiGp3CGl/Gsp5Twp5TVSyv8ppeyWUr4hpVwupbyi/Hcy01IAcZfvaoOOYssAA8k9ziSgdyoKube3t6Ovry906a7iXOSuXi1AVe4uRUZVPbodZBsY06dPBwAcPnzYqVJd5C6EGGR9xdmb0W0Zl90CDFTuNvvGRe7K+urs7HT2dUqdBE+umJQ2QFfuR48etba/XnfUPRzbnot+LNVnpe2KcylywC142tra0Nvbi1OnTln7cPA9P5yIexRyyMBVga5BK4QYdMzNRcauDVUVRyFtlzcbVIJxlHuYN2va8MnlcuzKXV/S2sj94MGDTlvm9OnT6O3tNbYrMPjkCtWWCatjPc41kWWzWTblPnXqVAC0Cc81aQM06yOfz1cu97gmdyEEjh49WrGWwsSHrqJd6amVLJUYTWMiisjytdsoq15T3/zCF76Aj3zkI6GfxcWwJ/eOjg50d3dbK7CtrQ1HjhwBYFf4wVMfrmUely0DDJwEXIrBRe5qg842MPRLW6Y6yWQylcmCk9xdg3vGjBkASuTusmWA0jI/LrmPHj0ahUJhwLnusLh8Po+WlhanjZbNZtHW1oZDhw45FWMmk0F7e7s1vWnTpgEokbtrz0Uvg6nuJk+ejKNHjwKgTQK29srn85gwYQKOHDli3RQOGzthccqqdLWrvg9hqpORI0eiqamJtKGq8kYld1u9BUWbKe4P/uAPcPPNN4d+FhfDmtzXrl2LCRMmYO3atdYKvPLKK3Hy5ElnZwmSe1hcY2MjmpubnaTd0tKCQqHgnASoy/Jx48Yhl8tVBpBpcCti3L17t9W+mTZt2oBNS1sndZG7Or7mIu1p06Yhm81iz5491qXvtGnTkM/nsW3bNqdyBy6paAq5m1Zk6mila6NMPddFFgAwe/Zs7Nu3z9pPstksJk6cOIAYbcr90KFD6O/vN7Zra2srent7cfLkSWs/mTx5ckXw2Npf1Z2tXVV6unIPe67qm3v37nVaUIDbgqSQuxCiMq5t7aAfwbS168iRI9Ha2lrpw6a2V6vP3bt3W+s3SQxrcl+4cCFyuRz+4z/+w1qBV155JQBg165d1gYJXlByxdlIW+9Uto48evRoNDc3Owd3NpvF9OnTsW/fPuugXbJkCQDgxRdftJLAjBkzsH//fgBuctc9dxeB2gZGoVDA7NmzsWPHDuuStlAoYOHChVi/fr3TXwYuqS1b3hQZuyZ31+kW4NKqxxU3Z86cAURmi3vttddI5L5//34r4c2fPx8AsGXLFmvclClTKqtAWz8J9mFbf9qzZ481btSoUZg6dSp27txptW/mzp0LANi2bZs1bxRy18vgaq8JEyYMsEdNcddeey3WrVtn7ZuzZs1CS0sLXn755ZTco2DUqFG49dZb8eSTT1pJe968eQCAl19+2Tq4FZH19/c7ScDluas4F7kDJYX36quvOtObNWuWk9znzZuH0aNHY8OGDdbBPWPGDBw/frxyucOl3G1kDAy83GWLmzt3LrZv3+6MW7ZsGTZt2oTu7m7rAAJKStA20Nra2tDf318Z4LZ2oJDx5Zdfjp07d5La9ciRIxX7wxR37bXXYuvWrRXCC3tuU1MT5syZU4kztesNN9wAoNTXbcQ4efJk9Pf3V45gmuImTpyI9vZ29Pb2AjCT+6JFi7Bz586KB24qq2p/m8KfN28e8vk8Nm/ejL6+PmNaynOXUjrJ3WV7AqVJdteuXc6+ftttt2H79u04dOiQMUYIgRtuuKFC7qZnJolhTe4AsGDBAuzduxc9PT3Gip4zZw5mzJiBn//851bSnjdvHrq7uytLaVNnmTNnzoDBTZ0ETHG33nornnvuuYoKMZVDLfNtHTmbzWLZsmVYv369ddDOnj0bALBjxw4rCUyfPh3Hjx+vvA7XlLepU6di3759TtK+6qqr8Nprr1XIwlQnb3nLW9DV1YVf//rX1qXvqFGj8Ktf/crartdeey0AOOOU0nbl7ZprrsGRI0cq7WWKe8c73gEA+M///E9r3IIFC3D27Fns378fmUzGWHcLFiyoKHIbkU2ePLlCjC6F/9JLL1nj5s6di3PnzlVWeaY4tWLcuHGjNW7BggXYunUrOjs7AYTXST6fx9VXX43Nmzdb+2ZbWxt6enoq+x+mOrnyyiuxffv2yruUTHHz58/H9u3bneP6+uuvBwBs377dStoLFy7E5s2brdyUJIY9uc+cORN9fX1ob2+3zqILFy6sqGNTg6gO/+qrrwIwN+7SpUtx9OhRbNmyBYC5s8yZMwc7duxAd3e3Ne7mm29GX18ftmzZYh3cs2bNQkdHB86ePWtMCygpy4MHD1oHrdrEeeqpp9Db22uMu+666yClxCuvvALATNo33HAD9u7dW1Fuprj58+ejp6cHu3fvtsapvYiOjg5jTCaTwdve9jY8//zzVnV04403IpfL4fnnn7e2/+WXX46enp7KXoRNaQPA5s2bK/kIw6JFi5DJZLBr1y5rnErv5Zdftrbr/PnzsXv3bnR1dRnbCyhNPtu2bbO2/5IlS5DP5/Hcc8+RVgIvvvgiADNpL168GACwYcMGAOa6W7hwIbq6uvDaa68hl8sZ2/a6666rrFJMaV1++eUASnarjdzf/va3o6urC7/85S8BmNth/vz5uHjxIvbu3WuN02+92kj7iiuuQFdXF44dO5aSexRMmTIFQGmjyTaLKo/RNrivvvpqAMDWrVsBmBv3tttuAwD87Gc/s8YtXboU58+fx/bt2wHY7RYA2LNnj3VwU+OmTp2KkydP4uzZs1bPde7cudi4cSP6+vqQz+dD49Tg3rRpEwAY4xYtWgSgpAQB+wACLk2gtmNuQOnYpK1dr732WuzcuRO9vb3GuBEjRmDhwoVOclckq4jMVQbVT0xtkc/nMXHiRKctc8UVVwAobb7Z2nX69OkoFos4dOiQldxnzpzpnNxHjBiBuXPnYseOHda4BQsWAEBlcrf532PGjKkQoyluzpw5AICdO3dayzp79mwcPXoUnZ2dzgnlueees5L7TTfdhGw2izVr1gCgt6trI//48eNW0lZ7Qq+//rpx3CSJYU/u6sjRhQsXrBWtzkR3d3cbG625uRnTpk1zDtpp06ahpaUFO3bssMapyWLbtm3WOP3EjI3IqOSunzawxc2ZMwf79++3KvepU6cin89jz549AMzkrvJ24MABAGbSVvsfrgGkXqp04sQJa7tOnToVvb29eOONN6x1t2DBAuzcudNK7osWLUJjYyOeffZZa96mTZuGQqHgVHjApUnKFqfK6lLkSsgcOHDAGjdt2jR0dHTgwoUL1rgZM2bgwIEDVutj1KhRaGtrw2uvvQbATNpAqd+5rA91pPPgwYOksh46dMh67Peaa67BL37xCyu5Nzc3Y9asWdi5cycAczvok6wtTr2grauri0Tub7zxBgqFgjEuKQx7ctcHj4vcgdJsaxuM8+fPdxIPUOqkishsJACUOqgtrrm5GblcztpBgUsECpgHD3Bp2ega3FOnTsXhw4etyk0Igba2Nhw+fBiAeXCrCUodrzOVddSoUZgyZUrlVFJDQ0NonPrqM9upJWDgiRlXex07dgxdXV3GuEKhUPGEbWVQF97U8UpbW1DIPZfLVb5315aWrgRt7ar6uituypQpznYFUDml5YpTosIWp/J27tw5Ulrt7e3WuAULFlROGtnqbsqUKU67bezYschkMpW+aeITdQkQMIsd4FJ7ueKSwpuO3G3KDShtviiVYoubPHlyZWDYdvPV2XRbnLpKbYtRZVAd3RWnYBsY6rSBzZYBULmMA5g76ciRI9Hc3FyZ8Kid3qRompqa0NTUBMDeriqtY8eOOdsLKPmkrtWRS7kBpTpWN4Gpyp0yCVDInZqWayWgv3PJNVmom7u2OCUqbPlraGhAS0uLNQYY+G1RLsFz4MAB9Pb2WuMmTZpUscdM7ZXJZAZ83R5FVDQ2NjrLAJj7eZIY9uQ+atSoSoezVaBOeLZGmzBhgvNIIlBqXHUyxBSXyWQqF1Rc6VEGtxCCFEcl9/HjxzvPLwM0cgdKA0ipI1tbKFXuSk9NeCZ1D1waZDa7DaC3/+TJkytfnuHqJ67TMgBNuetx1FWAa9KmxOmTBWUl4IqjKHfg0iRgi1FtT3lmsVjEuXPnYq+ggFKdqG9aovQnW9/Un5kq9wjQCc9W0dTBrXd4jri2tjbnBSDgUkdwnYflJHdq51MvQXLFTZo0qeK52shdDdxcLmctr4qzqSNqO0Rpf2odU0mFosptMcq+A4Y+uVNsQ456o7YDdz9Rcba+WSgUKl+WnZJ7RKiOQF0icSgtaqdSXrQrjkLaAI0ERo0aVZnoOEhAX27bOqleJxTl7lqqUgYQtb30MnC3f7WUuy5kqO1KJTyOSUCv47jKXV/dUfPGYWn5ijuboAQutUVqy0QERbnryzzOxnXF+So8F7lTSUA9t5qDW6/juHYLQCP3ESNGYMSIEQCqS8ZRlHtccqfGDRflbksrl8tVVC/3CorarpQ9PBdpqzpOlXtEUMhdvc0PKClbE7iVgK9idNkyFOWuP5fDlqHu+lM3kPSTMDaoMtjIHbg00GxlbW5uruSdSngcNl+pTpYCWIlXXmlyxKFyLd8VZytDY2NjZcIbyp67C5SyctsyukBRG/q29Lq6uowxelyq3COCQu4AKh1eKYIw1Fq5q9usJpQG0FLs23cH1q83x1EIL8oGHQe5qzj1JdgmqDK4JjwVZ5u0hRCVSUUdO7TlDRg40E3PdOWvvX0mgDUAPo93vCNjbLNSf1qK8+c/Zm1X1e9s7QXQiJGb3EukXZrIXnjB3NfVJKA2Lk2gnCDzPUHkitPbfPTo0ca4Up0sxaFD7yW1V6rcI4JK7opwbI2mN66aDMKgd3ibsvQj96Xo6voLa2fJ528CsAavvfY/sXw5jLEUW0YnOSoJcJB76dp4iQRsZVVksWvXH5AGkK1dAUDK0jNPnZpnjFGDFliJI0emO+JKsLX/0aNXAigAyKGnB3jmGVve1uDs2b+0tmtf3xIAK9HZeZ3xmQDQ0HALgJXlySUcpcmwVNZ9+9qMcXqdbNtmFkb790+Cmsh++7eFsQyqXU+e/LC1XSkTVKFQqEzq1EmAstkP2Md/V9cNKLXXJ6ztlXruMaFus7m+i1B1EptyLzVCqSOfOHGFMU4pLWAljh6dYYzTB4ZtWV4ajGtw4cJKa2cpFm8C0AAgayWL11+fBWAlenoWGZ9ZUpylvB04MMkYp5f1179uMcZRPffe3sVQJGAra2fn9QDW4OjR/0WayGzKff164MSJxwF8Hg899PvGtPbunVjJ2z33zHc8s1QnBw9OMT53xYpGNDZmkM1KFArALbeEx508uQBCNMI2CaxfDzz55McBfB5r1/61MW/r1wN79jwM4PN4/PG7jHEbNohKWVet+m1jXKl8pbiPfGSOMW7jxiYI0WAtAwB0dFwOSvuX+tNS7Nxpbi8AaGz8LbjGoW6PlcZGOEoXD0txpfoxlWE+MpkmuMah4hp1z6KqUN/YUss/ixYtknHw4osvSgDyc5/7nDXuT//0TyUA+dJLLxlj1q2TErgggV7Z0NAn160Lj1u7trcS19RUNMY9+uh2LU4a4x58UMpMpigBKbPZ0r9N+WtqKsWY0lu3Tsp8vsdZhnXrpBTiojNuzZqLWp30G+M2bNgggaUSWCl/8Yuu8CCPst5/f78Uot8Z9z/+x99LYKV873u/bn2mEH3ltIrGtEpx7mf+13+dJbW/lKV6fvBBc9urGFe7Dqw3ahm44tx1RymDlFKuXHlGAr3OOr7jjn8o13Gfta+rdigUeo3PfP75YiWusdHch//930+Qxiu1rG+88Yb87Gc/K8+cORMeEBMANkkDr9ac2CUDuUsp5csvvyzPnTtnjenq6pIvvviiNcaHBDIZNwl8/vN9pI5M7Swq1kYWDz5YyrvrmQPLYC/rpTKY455++jyJ8KhlpU5khUKvBHplPt/jSKsos9kiy6B94IEiqV194GpX/3rjKStn37yUnrstPvzhvc469pmgAOq4do8dalmTxpuC3LnATQLU9FQsR2fhHrTUMlBVr0qTUla/icw8aDmfqWKo7coJzjIkEUcFtY4bG+117N+H+SayoQAbuQtZvj1ZSyxevFiqV8oOBaxfX/LQbrkFWLasenGcqEUZ1q8Hli8HenqAQgFYsyb58tbimfqzq92ubzZQ+91wH69RIYT4lZRycehnKbmn4MRQnshSpKg3pOSeIkWKFHUIG7nXxVHIFClSpEgxECm5p0iRIkUdIiX3FClSpKhDpOSeIkWKFHWIlNxTpEiRog6RknuKFClS1CGGxFFIIUQHgAMxkhgP4ARTdjiR5ssfQzVvab78MVTzVk/5miGlbA37YEiQe1wIITaZznrWEmm+/DFU85bmyx9DNW9vlnyltkyKFClS1CFSck+RIkWKOkS9kPsjtc6AAWm+/DFU85bmyx9DNW9vinzVheeeIkWKFCkGol6Ue4oUKVKk0JCSe4oUKVLUIYYNuQshVgghdgohdgshVoZ8LoQQXy1/vkUIsbBK+ZomhPiFEGK7EGKbEOKekJhbhBBnhBCvlP98pkp52y+E2Fp+5qB3KteizoQQc7V6eEUIcVYI8eeBmKrVlxDiUSHEcSHEq9rvxgohnhZC7Cr/Pcbwf619MoF8/Y0QYke5rX4khLjM8H+t7Z5Q3v6PEOKI1ma/Y/i/1a6z72t52i+EeMXwfxOrMxNHJN7PTF/RNJT+AMgC2ANgNoACgM0Arg7E/A6AJwEIlL6+fGOV8jYJwMLyz6MBvBaSt1sA/KQG9bYfwHjL5zWps0C7vo7SRYya1BeAmwAsBPCq9rsvAVhZ/nklgC8a8m7tkwnk678ByJV//mJYvijtnlDe/g+AvyS0d1XrLPD53wL4TLXrzMQRSfez4aLcbwSwW0q5V0rZA+BxAO8KxLwLwLdlCRsAXCaEmJR0xqSU7VLKl8o/nwOwHcCUpJ/LhJrUmYblAPZIKePcTo4FKeVaACcDv34XgMfKPz8G4PaQ/0rpk6z5klI+JaXsK/9zA4CpXM/zgaHOKKh6nSkIIQSAPwLwPa7nUWHhiET72XAh9ykADmn/PozBBEqJSRRCiJkAbgCwMeTjZUKIzUKIJ4UQ86uUJQngKSHEr4QQHw75vNZ1dgfMg60W9aUwUUrZDpQGJoAJITG1rrsPoLTqCoOr3ZPC/y5bRo8aLIZa1tnbAByTUu4yfF6VOgtwRKL9bLiQuwj5XfAMJyUmMQghRgH4FwB/LqU8G/j4JZSsh+sAfA3Av1UpW2+VUi4E8E4Adwshbgp8XrM6E0IUAPwegP8v5ONa1ZcPall39wHoA/AdQ4ir3ZPANwDMAXA9gHaULJAgajlG3wO7ak+8zhwcYfxvIb8j1dlwIffDAKZp/54K4GiEmEQghMij1GjfkVL+a/BzKeVZKeX58s//CSAvhBifdL6klEfLfx8H8COUlng6alZnKA2il6SUx4If1Kq+NBxT9lT57+MhMTWpOyHEnQD+O4A/lmVTNghCu7NDSnlMStkvpSwC+CfDM2tVZzkAvw/g+6aYpOvMwBGJ9rPhQu4vArhCCDGrrPjuAPBEIOYJAO8rnwBZCuCMWvIkibKX900A26WUf2eIaSvHQQhxI0r1/kbC+RophBitfkZpM+7VQFhN6qwMo5KqRX0F8ASAO8s/3wngxyExlD7JCiHECgCfBPB7UsqLhhhKuyeRN32v5t2GZ1a9zsr4bQA7pJSHwz5Mus4sHJFsP0tidzihHeffQWmXeQ+A+8q/+yiAj5Z/FgC+Xv58K4DFVcrXb6K0TNoC4JXyn98J5O1/A9iG0k73BgC/UYV8zS4/b3P52UOpzkagRNYt2u9qUl8oTTDtAHpRUkkfBDAOwBoAu8p/jy3HTgbwn7Y+mXC+dqPkv6p+9o/BfJnavQp5W13uQ1tQIp9JQ6HOyr//lupbWmzV6szCEYn2s/T1AylSpEhRhxgutkyKFClSpPBASu4pUqRIUYdIyT1FihQp6hApuadIkSJFHSIl9xQpUqSoQ6TkniJFihR1iJTcU6RIkaIO8f8DVtdTE6SUErAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pks = fidp[\"pks\"]\n",
    "ons = fidp[\"ons\"]\n",
    "t = np.arange(0,(len(abp)/fs),1.0/fs)\n",
    "plt.plot(t, abp, color = 'black')\n",
    "plt.plot(t[pks], abp[pks], \".\", color = 'red')\n",
    "plt.plot(t[ons], abp[ons], \".\", color = 'blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate BP values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Extract systolic and diastolic BPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "sbp = np.median(abp[fidp['pks']])\n",
    "dbp = np.median(abp[fidp['ons']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Extract mean BPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "ons = fidp['ons']\n",
    "off = fidp['off']\n",
    "mbps = np.zeros(len(ons))\n",
    "for beat_no in range(0,len(ons)):\n",
    "    mbps[beat_no] = np.mean(abp[ons[beat_no]:off[beat_no]])\n",
    "mbp = np.median(mbps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Print results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Systolic blood pressure  (SBP): 161.9 mmHg\n",
      "Diastolic blood pressure (DBP): 67.2 mmHg\n",
      "Mean blood pressure      (MBP): 102.5 mmHg\n"
     ]
    }
   ],
   "source": [
    "print('Systolic blood pressure  (SBP): {:.1f} mmHg'.format(sbp))\n",
    "print('Diastolic blood pressure (DBP): {:.1f} mmHg'.format(dbp))\n",
    "print('Mean blood pressure      (MBP): {:.1f} mmHg'.format(mbp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\"> <b>Question:</b>\n",
    "    Often MBP is approximated as: MBP = ((2/3)*DBP)+((1/3)*SBP).\n",
    "    Would it have made much difference if we had used this approximation?\n",
    "    What is the basis for this approximation?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Beat Detection Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following functions are required for this tutorial. Run the cell below and then return to the top of the page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import scipy.signal as sp\n",
    "import numpy as np\n",
    "\n",
    "def pulse_detect(x,fs,w,alg):\n",
    "    \"\"\"\n",
    "    Description: Pulse detection and correction from pulsatile signals\n",
    "    Inputs:  x, array with pulsatile signal [user defined units]\n",
    "             fs, sampling rate of signal [Hz]\n",
    "             w, window length for analysis [s]\n",
    "             alg, string with the name of the algorithm to apply ['heartpy','d2max','upslopes','delineator']\n",
    "    Outputs: ibis, location of cardiac cycles as detected by the selected algorithm [number of samples]\n",
    "\n",
    "    Algorithms:       1: HeartPy (van Gent et al, 2019, DOI: 10.1016/j.trf.2019.09.015)\n",
    "                      2: 2nd derivative maxima (Elgendi et al, 2013, DOI: 10.1371/journal.pone.0076585)\n",
    "                      3: Systolic upslopes (Arguello Prada and Serna Maldonado, 2018,\n",
    "                         DOI: 10.1080/03091902.2019.1572237)\n",
    "                      4: Delineator (Li et al, 2010, DOI: 10.1109/TBME.2005.855725)\n",
    "    Fiducial points:  1: Systolic peak (pks)\n",
    "                      2: Onset, as the minimum before the systolic peak (ons)\n",
    "                      3: Onset, using the tangent intersection method (ti)\n",
    "                      4: Diastolic peak (dpk)\n",
    "                      5: Maximum slope (m1d)\n",
    "                      6: a point from second derivative PPG (a2d)\n",
    "                      7: b point from second derivative PPG (b2d)\n",
    "                      8: c point from second derivative PPG (c2d)\n",
    "                      9: d point from second derivative PPG (d2d)\n",
    "                      10: e point from second derivative PPG (e2d)\n",
    "                      11: p1 from the third derivative PPG (p1)\n",
    "                      12: p2 from the third derivative PPG (p2)\n",
    "\n",
    "    Libraries: NumPy (as np), SciPy (Signal, as sp), Matplotlib (PyPlot, as plt)\n",
    "\n",
    "    Version: 1.0 - June 2022\n",
    "\n",
    "    Developed by: Elisa Mejía-Mejía\n",
    "                   City, University of London\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Check selected algorithm\n",
    "    pos_alg = ['heartpy','d2max','upslopes','delineator']\n",
    "    if not(alg in pos_alg):\n",
    "        print('Unknown algorithm determined. Using D2max as default')\n",
    "        alg = 'd2max'\n",
    "\n",
    "    # Pre-processing of signal\n",
    "    x_d = sp.detrend(x)\n",
    "    sos = sp.butter(10, [0.5, 10], btype = 'bp', analog = False, output = 'sos', fs = fs)\n",
    "    x_f = sp.sosfiltfilt(sos, x_d)\n",
    "\n",
    "    # Peak detection in windows of length w\n",
    "    n_int = np.floor(len(x_f)/(w*fs))\n",
    "    for i in range(int(n_int)):\n",
    "        start = i*fs*w\n",
    "        stop = (i + 1)*fs*w - 1\n",
    "        # print('Start: ' + str(start) + ', stop: ' + str(stop) + ', fs: ' + str(fs))\n",
    "        aux = x_f[range(start,stop)]\n",
    "        if alg == 'heartpy':\n",
    "            locs = heartpy(aux,fs,40,180,5)\n",
    "        elif alg == 'd2max':\n",
    "            locs = d2max(aux,fs)\n",
    "        elif alg == 'upslopes':\n",
    "            locs = upslopes(aux)\n",
    "        elif alg == 'delineator':\n",
    "            locs = delineator(aux,fs)\n",
    "        locs = locs + start\n",
    "        if i == 0:\n",
    "            ibis = locs\n",
    "        else:\n",
    "            ibis = np.append(ibis,locs)\n",
    "    if n_int*fs*w != len(x_f):\n",
    "        start = stop + 1\n",
    "        stop = len(x_f)\n",
    "        aux = x_f[range(start,stop)]\n",
    "        if len(aux) > 20:\n",
    "            if alg == 'heartpy':\n",
    "                locs = heartpy(aux,fs,40,180,5)\n",
    "            elif alg == 'd2max':\n",
    "                locs = d2max(aux,fs)\n",
    "            elif alg == 'upslopes':\n",
    "                locs = upslopes(aux)\n",
    "            elif alg == 'delineator':\n",
    "                locs = delineator(aux,fs)\n",
    "            locs = locs + start\n",
    "            ibis = np.append(ibis,locs)\n",
    "    ind, = np.where(ibis <= len(x_f))\n",
    "    ibis = ibis[ind]\n",
    "\n",
    "    ibis = peak_correction(x,ibis,fs,20,5,[0.5, 1.5])\n",
    "\n",
    "    #fig = plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.plot(x_d)\n",
    "    #plt.plot(x_f)\n",
    "    #plt.scatter(ibis,x_f[ibis],marker = 'o',color = 'red')\n",
    "    #plt.scatter(ibis,x[ibis],marker = 'o',color = 'red')\n",
    "\n",
    "    return ibis\n",
    "\n",
    "def peak_correction(x,locs,fs,t,stride,th_len):\n",
    "    \"\"\"\n",
    "    Correction of peaks detected from pulsatile signals\n",
    "\n",
    "    Inputs:   x, pulsatile signal [user defined units]\n",
    "              locs, location of the detected interbeat intervals [number of samples]\n",
    "              fs, sampling rate [Hz]\n",
    "              t, duration of intervals for the correction [s]\n",
    "              stride, stride between consecutive intervals for the correction [s]\n",
    "              th_len, array with the percentage of lower and higher thresholds for comparing the duration of IBIs\n",
    "              [proportions]\n",
    "    Outputs:  ibis, array with the corrected points related to the start of the inter-beat intervals [number of samples]\n",
    "\n",
    "    Developed by:  Elisa Mejía Mejía\n",
    "                   City, University of London\n",
    "    Version:       1.0 -   June, 2022\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    #fig = plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.scatter(locs,x[locs],marker = 'o',color = 'red', label = 'Original')\n",
    "    #plt.title('Peak correction')\n",
    "\n",
    "    # Correction of long and short IBIs\n",
    "    len_window = np.round(t*fs)\n",
    "    #print('Window length: ' + str(len_window))\n",
    "    first_i = 0\n",
    "    second_i = len_window - 1\n",
    "    while second_i < len(x):\n",
    "        ind1, = np.where(locs >= first_i)\n",
    "        ind2, = np.where(locs <= second_i)\n",
    "        ind = np.intersect1d(ind1, ind2)\n",
    "\n",
    "        win = locs[ind]\n",
    "        dif = np.diff(win)\n",
    "        #print('Indices: ' + str(ind) + ', locs: ' + str(locs[ind]) + ', dif: ' + str(dif))\n",
    "\n",
    "        th_dif = np.zeros(2)\n",
    "        th_dif[0] = th_len[0]*np.median(dif)\n",
    "        th_dif[1] = th_len[1]*np.median(dif)\n",
    "\n",
    "        th_amp = np.zeros(2)\n",
    "        th_amp[0] = 0.75*np.median(x[win])\n",
    "        th_amp[1] = 1.25*np.median(x[win])\n",
    "        #print('Length thresholds: ' + str(th_dif) + ', amplitude thresholds: ' + str(th_amp))\n",
    "\n",
    "        j = 0\n",
    "        while j < len(dif):\n",
    "            if dif[j] <= th_dif[0]:\n",
    "                if j == 0:\n",
    "                    opt = np.append(win[j], win[j + 1])\n",
    "                else:\n",
    "                    opt = np.append(win[j], win[j + 1]) - win[j - 1]\n",
    "                print('Optional: ' + str(opt))\n",
    "                dif_abs = np.abs(opt - np.median(dif))\n",
    "                min_val = np.min(dif_abs)\n",
    "                ind_min, = np.where(dif_abs == min_val)\n",
    "                print('Minimum: ' + str(min_val) + ', index: ' + str(ind_min))\n",
    "                if ind_min == 0:\n",
    "                    print('Original window: ' + str(win), end = '')\n",
    "                    win = np.delete(win, win[j + 1])\n",
    "                    print(', modified window: ' + str(win))\n",
    "                else:\n",
    "                    print('Original window: ' + str(win), end = '')\n",
    "                    win = np.delete(win, win[j])\n",
    "                    print(', modified window: ' + str(win))\n",
    "                dif = np.diff(win)\n",
    "            elif dif[j] >= th_dif[1]:\n",
    "                aux_x = x[win[j]:win[j + 1]]\n",
    "                locs_pks, _ = sp.find_peaks(aux_x)\n",
    "                #fig = plt.figure()\n",
    "                #plt.plot(aux_x)\n",
    "                #plt.scatter(locs_pks,aux_x[locs_pks],marker = 'o',color = 'red')\n",
    "\n",
    "                locs_pks = locs_pks + win[j]\n",
    "                ind1, = np.where(x[locs_pks] >= th_amp[0])\n",
    "                ind2, = np.where(x[locs_pks] <= th_amp[1])\n",
    "                ind = np.intersect1d(ind1, ind2)\n",
    "                locs_pks = locs_pks[ind]\n",
    "                #print('Locations: ' + str(locs_pks))\n",
    "\n",
    "                if len(locs_pks) != 0:\n",
    "                    opt = locs_pks - win[j]\n",
    "\n",
    "                    dif_abs = np.abs(opt - np.median(dif))\n",
    "                    min_val = np.min(dif_abs)\n",
    "                    ind_min, = np.where(dif_abs == min_val)\n",
    "\n",
    "                    win = np.append(win, locs_pks[ind_min])\n",
    "                    win = np.sort(win)\n",
    "                    dif = np.diff(win)\n",
    "                    j = j + 1\n",
    "                else:\n",
    "                    opt = np.round(win[j] + np.median(dif))\n",
    "                    if opt < win[j + 1]:\n",
    "                        win = np.append(win, locs_pks[ind_min])\n",
    "                        win = np.sort(win)\n",
    "                        dif = np.diff(win)\n",
    "                        j = j + 1\n",
    "                    else:\n",
    "                        j = j + 1\n",
    "            else:\n",
    "                j = j + 1\n",
    "\n",
    "        locs = np.append(win, locs)\n",
    "        locs = np.sort(locs)\n",
    "\n",
    "        first_i = first_i + stride*fs - 1\n",
    "        second_i = second_i + stride*fs - 1\n",
    "\n",
    "    dif = np.diff(locs)\n",
    "    dif = np.append(0, dif)\n",
    "    ind, = np.where(dif != 0)\n",
    "    locs = locs[ind]\n",
    "\n",
    "    #plt.scatter(locs,x[locs],marker = 'o',color = 'green', label = 'After length correction')\n",
    "\n",
    "    # Correction of points that are not peaks\n",
    "    i = 0\n",
    "    pre_loc = 0\n",
    "    while i < len(locs):\n",
    "        if locs[i] == 0:\n",
    "            locs = np.delete(locs, locs[i])\n",
    "        elif locs[i] == len(x):\n",
    "            locs = np.delete(locs, locs[i])\n",
    "        else:\n",
    "            #print('Previous: ' + str(x[locs[i] - 1]) + ', actual: ' + str(x[locs[i]]) + ', next: ' + str(x[locs[i] + 1]))\n",
    "            cond = (x[locs[i]] >= x[locs[i] - 1]) and (x[locs[i]] >= x[locs[i] + 1])\n",
    "            #print('Condition: ' + str(cond))\n",
    "            if cond:\n",
    "                i = i + 1\n",
    "            else:\n",
    "                if locs[i] == pre_loc:\n",
    "                    i = i + 1\n",
    "                else:\n",
    "                    if i == 0:\n",
    "                        aux = x[0:locs[i + 1] - 1]\n",
    "                        aux_loc = locs[i] - 1\n",
    "                        aux_start = 0\n",
    "                    elif i == len(locs) - 1:\n",
    "                        aux = x[locs[i - 1]:len(x) - 1]\n",
    "                        aux_loc = locs[i] - locs[i - 1]\n",
    "                        aux_start = locs[i - 1]\n",
    "                    else:\n",
    "                        aux = x[locs[i - 1]:locs[i + 1]]\n",
    "                        aux_loc = locs[i] - locs[i - 1]\n",
    "                        aux_start = locs[i - 1]\n",
    "                    #print('i ' + str(i) + ' out of ' + str(len(locs)) + ', aux length: ' + str(len(aux)) +\n",
    "                    #      ', location: ' + str(aux_loc))\n",
    "                    #print('Locs i - 1: ' + str(locs[i - 1]) + ', locs i: ' + str(locs[i]) + ', locs i + 1: ' + str(locs[i + 1]))\n",
    "\n",
    "                    pre = find_closest_peak(aux, aux_loc, 'backward')\n",
    "                    pos = find_closest_peak(aux, aux_loc, 'forward')\n",
    "                    #print('Previous: ' + str(pre) + ', next: ' + str(pos) + ', actual: ' + str(aux_loc))\n",
    "\n",
    "                    ibi_pre = np.append(pre - 1, len(aux) - pre)\n",
    "                    ibi_pos = np.append(pos - 1, len(aux) - pos)\n",
    "                    ibi_act = np.append(aux_loc - 1, len(aux) - aux_loc)\n",
    "                    #print('Previous IBIs: ' + str(ibi_pre) + ', next IBIs: ' + str(ibi_pos) +\n",
    "                    #      ', actual IBIs: ' + str(ibi_act))\n",
    "\n",
    "                    dif_pre = np.abs(ibi_pre - np.mean(np.diff(locs)))\n",
    "                    dif_pos = np.abs(ibi_pos - np.mean(np.diff(locs)))\n",
    "                    dif_act = np.abs(ibi_act - np.mean(np.diff(locs)))\n",
    "                    #print('Previous DIF: ' + str(dif_pre) + ', next DIF: ' + str(dif_pos) +\n",
    "                    #      ', actual DIF: ' + str(dif_act))\n",
    "\n",
    "                    avgs = [np.mean(dif_pre), np.mean(dif_pos), np.mean(dif_act)]\n",
    "                    min_avg = np.min(avgs)\n",
    "                    ind, = np.where(min_avg == avgs)\n",
    "                    #print('Averages: ' + str(avgs) + ', min index: ' + str(ind))\n",
    "                    if len(ind) != 0:\n",
    "                        ind = ind[0]\n",
    "\n",
    "                    if ind == 0:\n",
    "                        locs[i] = pre + aux_start - 1\n",
    "                    elif ind == 1:\n",
    "                        locs[i] = pos + aux_start - 1\n",
    "                    elif ind == 2:\n",
    "                        locs[i] = aux_loc + aux_start - 1\n",
    "                    i = i + 1\n",
    "\n",
    "    #plt.scatter(locs,x[locs],marker = 'o',color = 'yellow', label = 'After not-peak correction')\n",
    "\n",
    "    # Correction of peaks according to amplitude\n",
    "    len_window = np.round(t*fs)\n",
    "    #print('Window length: ' + str(len_window))\n",
    "    keep = np.empty(0)\n",
    "    first_i = 0\n",
    "    second_i = len_window - 1\n",
    "    while second_i < len(x):\n",
    "        ind1, = np.where(locs >= first_i)\n",
    "        ind2, = np.where(locs <= second_i)\n",
    "        ind = np.intersect1d(ind1, ind2)\n",
    "        win = locs[ind]\n",
    "        if np.median(x[win]) > 0:\n",
    "            th_amp_low = 0.5*np.median(x[win])\n",
    "            th_amp_high = 3*np.median(x[win])\n",
    "        else:\n",
    "            th_amp_low = -3*np.median(x[win])\n",
    "            th_amp_high = 1.5*np.median(x[win])\n",
    "        ind1, = np.where(x[win] >= th_amp_low)\n",
    "        ind2, = np.where(x[win] <= th_amp_high)\n",
    "        aux_keep = np.intersect1d(ind1,ind2)\n",
    "        keep = np.append(keep, aux_keep)\n",
    "\n",
    "        first_i = second_i + 1\n",
    "        second_i = second_i + stride*fs - 1\n",
    "\n",
    "    if len(keep) != 0:\n",
    "        keep = np.unique(keep)\n",
    "        locs = locs[keep.astype(int)]\n",
    "\n",
    "    #plt.scatter(locs,x[locs],marker = 'o',color = 'purple', label = 'After amplitude correction')\n",
    "    #plt.legend()\n",
    "\n",
    "    return locs\n",
    "\n",
    "def find_closest_peak(x, loc, dir_search):\n",
    "    \"\"\"\n",
    "    Finds the closest peak to the initial location in x\n",
    "\n",
    "    Inputs:   x, signal of interest [user defined units]\n",
    "              loc, initial location [number of samples]\n",
    "              dir_search, direction of search ['backward','forward']\n",
    "    Outputs:  pos, location of the first peak detected in specified direction [number of samples]\n",
    "\n",
    "    Developed by:  Elisa Mejía Mejía\n",
    "                   City, University of London\n",
    "    Version:       1.0 -   June, 2022\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    pos = -1\n",
    "    if dir_search == 'backward':\n",
    "        i = loc - 2\n",
    "        while i > 0:\n",
    "            if (x[i] > x[i - 1]) and (x[i] > x[i + 1]):\n",
    "                pos = i\n",
    "                i = 0\n",
    "            else:\n",
    "                i = i - 1\n",
    "        if pos == -1:\n",
    "            pos = loc\n",
    "    elif dir_search == 'forward':\n",
    "        i = loc + 1\n",
    "        while i < len(x) - 1:\n",
    "            if (x[i] > x[i - 1]) and (x[i] > x[i + 1]):\n",
    "                pos = i\n",
    "                i = len(x)\n",
    "            else:\n",
    "                i = i + 1\n",
    "        if pos == -1:\n",
    "            pos = loc\n",
    "\n",
    "    return pos\n",
    "\n",
    "def seek_local(x, start, end):\n",
    "    val_min = x[start]\n",
    "    val_max = x[start]\n",
    "\n",
    "    ind_min = start\n",
    "    ind_max = start\n",
    "\n",
    "    for j in range(start, end):\n",
    "        if x[j] > val_max:\n",
    "            val_max = x[j]\n",
    "            ind_max = j\n",
    "        elif x[j] < val_min:\n",
    "            val_min = x[j]\n",
    "            ind_min = j\n",
    "\n",
    "    return val_min, ind_min, val_max, ind_max\n",
    "\n",
    "def heartpy(x, fs, min_ihr, max_ihr, w):\n",
    "    \"\"\"\n",
    "    Detects inter-beat intervals using HeartPy\n",
    "    Citation: van Gent P, Farah H, van Nes N, van Arem B (2019) Heartpy: A novel heart rate algorithm\n",
    "              for the analysis of noisy signals. Transp Res Part F, vol. 66, pp. 368-378. DOI: 10.1016/j.trf.2019.09.015\n",
    "\n",
    "    Inputs:   x, pulsatile signal [user defined units]\n",
    "              fs, sampling rate [Hz]\n",
    "              min_ihr, minimum value of instantaneous heart rate to be accepted [bpm]\n",
    "              max_ihr, maximum value of instantaneous heart rate to be accepted [bpm]\n",
    "              w, length of segments for correction of peaks [s]\n",
    "    Outputs:  ibis, position of the starting points of inter-beat intervals [number of samples]\n",
    "\n",
    "    Developed by:  Elisa Mejía Mejía\n",
    "                   City, University of London\n",
    "    Version:       1.0 -   June, 2022\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Identification of peaks\n",
    "    is_roi = 0\n",
    "    n_rois = 0\n",
    "    pos_pks = np.empty(0).astype(int)\n",
    "    locs = np.empty(0).astype(int)\n",
    "\n",
    "    len_ma = int(np.round(0.75*fs))\n",
    "    #print(len_ma)\n",
    "    sig = np.append(x[0]*np.ones(len_ma), x)\n",
    "    sig = np.append(sig, x[-1]*np.ones(len_ma))\n",
    "\n",
    "    i = len_ma\n",
    "    while i < len(sig) - len_ma:\n",
    "        ma = np.mean(sig[i - len_ma:i + len_ma - 1])\n",
    "        #print(len(sig[i - len_ma:i + len_ma - 1]),ma)\n",
    "\n",
    "        # If it is the beginning of a new ROI:\n",
    "        if is_roi == 0 and sig[i] >= ma:\n",
    "            is_roi = 1\n",
    "            n_rois = n_rois + 1\n",
    "            #print('New ROI ---' + str(n_rois) + ' @ ' + str(i))\n",
    "            # If it is a peak:\n",
    "            if sig[i] >= sig[i - 1] and sig[i] >= sig[i + 1]:\n",
    "                pos_pks = np.append(pos_pks, int(i))\n",
    "                #print('Possible peaks: ' + str(pos_pks))\n",
    "\n",
    "        # If it is part of a ROI which is not over:\n",
    "        elif is_roi == 1 and sig[i] > ma:\n",
    "            #print('Actual ROI ---' + str(n_rois) + ' @ ' + str(i))\n",
    "            # If it is a peak:\n",
    "            if sig[i] >= sig[i - 1] and sig[i] >= sig[i + 1]:\n",
    "                pos_pks = np.append(pos_pks, int(i))\n",
    "                #print('Possible peaks: ' + str(pos_pks))\n",
    "\n",
    "        # If the ROI is over or the end of the signal has been reached:\n",
    "        elif is_roi == 1 and (sig[i] < ma or i == (len(sig) - len_ma)):\n",
    "            #print('End of ROI ---' + str(n_rois) + ' @ ' + str(i) + '. Pos pks: ' + str(pos_pks))\n",
    "            is_roi = 0 # Lowers flag\n",
    "\n",
    "            # If it is the end of the first ROI:\n",
    "            if n_rois == 1:\n",
    "                # If at least one peak has been found:\n",
    "                if len(pos_pks) != 0:\n",
    "                    # Determines the location of the maximum peak:\n",
    "                    max_pk = np.max(sig[pos_pks])\n",
    "                    ind, = np.where(max_pk == np.max(sig[pos_pks]))\n",
    "                    #print('First ROI: (1) Max Peak: ' + str(max_pk) + ', amplitudes: ' + str(sig[pos_pks]) +\n",
    "                    #      ', index: ' + str(int(ind)), ', pk_ind: ' + str(pos_pks[ind]))\n",
    "                    # The maximum peak is added to the list:\n",
    "                    locs = np.append(locs, pos_pks[ind])\n",
    "                    #print('Locations: ' + str(locs))\n",
    "                # If no peak was found:\n",
    "                else:\n",
    "                    # Counter for ROIs is reset to previous value:\n",
    "                    n_rois = n_rois - 1\n",
    "\n",
    "            # If it is the end of the second ROI:\n",
    "            elif n_rois == 2:\n",
    "                # If at least one peak has been found:\n",
    "                if len(pos_pks) != 0:\n",
    "                    # Measures instantantaneous HR of found peaks with respect to the previous peak:\n",
    "                    ihr = 60/((pos_pks - locs[-1])/fs)\n",
    "                    good_ihr, = np.where(ihr <= max_ihr and ihr >= min_ihr)\n",
    "                    #print('Second ROI IHR check: (1) IHR: ' + str(ihr) + ', valid peaks: ' + str(good_ihr) +\n",
    "                    #      ', pos_pks before: ' + str(pos_pks) + ', pos_pks after: ' + str(pos_pks[good_ihr]))\n",
    "                    pos_pks = pos_pks[good_ihr].astype(int)\n",
    "\n",
    "                    # If at least one peak is between HR limits:\n",
    "                    if len(pos_pks) != 0:\n",
    "                        # Determines the location of the maximum peak:\n",
    "                        max_pk = np.max(sig[pos_pks])\n",
    "                        ind, = np.where(max_pk == np.max(sig[pos_pks]))\n",
    "                        #print('Second ROI: (1) Max Peak: ' + str(max_pk) + ', amplitudes: ' + str(sig[pos_pks]) +\n",
    "                        #  ', index: ' + str(int(ind)), ', pk_ind: ' + str(pos_pks[ind]))\n",
    "                        # The maximum peak is added to the list:\n",
    "                        locs = np.append(locs, pos_pks[ind])\n",
    "                        #print('Locations: ' + str(locs))\n",
    "                # If no peak was found:\n",
    "                else:\n",
    "                    # Counter for ROIs is reset to previous value:\n",
    "                    n_rois = n_rois - 1\n",
    "\n",
    "            # If it is the end of the any further ROI:\n",
    "            else:\n",
    "                # If at least one peak has been found:\n",
    "                if len(pos_pks) != 0:\n",
    "                    # Measures instantantaneous HR of found peaks with respect to the previous peak:\n",
    "                    ihr = 60/((pos_pks - locs[-1])/fs)\n",
    "                    good_ihr, = np.where(ihr <= max_ihr and ihr >= min_ihr)\n",
    "                    #print('Third ROI IHR check: (1) IHR: ' + str(ihr) + ', valid peaks: ' + str(good_ihr) +\n",
    "                    #      ', pos_pks before: ' + str(pos_pks) + ', pos_pks after: ' + str(pos_pks[good_ihr]))\n",
    "                    pos_pks = pos_pks[good_ihr].astype(int)\n",
    "\n",
    "                    # If at least one peak is between HR limits:\n",
    "                    if len(pos_pks) != 0:\n",
    "                        # Calculates SDNN with the possible peaks on the ROI:\n",
    "                        sdnn = np.zeros(len(pos_pks))\n",
    "                        for j in range(len(pos_pks)):\n",
    "                            sdnn[j] = np.std(np.append(locs/fs, pos_pks[j]/fs))\n",
    "                        # Determines the new peak as that one with the lowest SDNN:\n",
    "                        min_pk = np.min(sdnn)\n",
    "                        ind, = np.where(min_pk == np.min(sdnn))\n",
    "                        #print('Third ROI: (1) Min SDNN Peak: ' + str(min_pk) + ', amplitudes: ' + str(sig[pos_pks]) +\n",
    "                        #  ', index: ' + str(int(ind)), ', pk_ind: ' + str(pos_pks[ind]))\n",
    "                        locs = np.append(locs, pos_pks[ind])\n",
    "                        #print('Locations: ' + str(locs))\n",
    "                # If no peak was found:\n",
    "                else:\n",
    "                    # Counter for ROIs is reset to previous value:\n",
    "                    n_rois = n_rois - 1\n",
    "\n",
    "            # Resets possible peaks for next ROI:\n",
    "            pos_pks = np.empty(0)\n",
    "\n",
    "        i = i + 1;\n",
    "\n",
    "    locs = locs - len_ma\n",
    "\n",
    "    # Correction of peaks\n",
    "    c_locs = np.empty(0)\n",
    "    n_int = np.floor(len(x)/(w*fs))\n",
    "    for i in range(int(n_int)):\n",
    "        ind1, = np.where(locs >= i*w*fs)\n",
    "        #print('Locs >= ' + str((i)*w*fs) + ': ' + str(locs[ind1]))\n",
    "        ind2, = np.where(locs < (i + 1)*w*fs)\n",
    "        #print('Locs < ' + str((i + 1)*w*fs) + ': ' + str(locs[ind2]))\n",
    "        ind = np.intersect1d(ind1, ind2)\n",
    "        #print('Larger and lower than locs: ' + str(locs[ind]))\n",
    "        int_locs = locs[ind]\n",
    "\n",
    "        if i == 0:\n",
    "            aux_ibis = np.diff(int_locs)\n",
    "        else:\n",
    "            ind, = np.where(locs >= i*w*fs)\n",
    "            last = locs[ind[0] - 1]\n",
    "            aux_ibis = np.diff(np.append(last, int_locs))\n",
    "        avg_ibis = np.mean(aux_ibis)\n",
    "        th = np.append((avg_ibis - 0.3*avg_ibis), (avg_ibis + 0.3*avg_ibis))\n",
    "        ind1, = np.where(aux_ibis > th[0])\n",
    "        #print('Ind1: ' + str(ind1))\n",
    "        ind2, = np.where(aux_ibis < th[1])\n",
    "        #print('Ind2: ' + str(ind2))\n",
    "        ind = np.intersect1d(ind1, ind2)\n",
    "        #print('Ind: ' + str(ind))\n",
    "\n",
    "        c_locs = np.append(c_locs, int_locs[ind]).astype(int)\n",
    "        print(c_locs)\n",
    "\n",
    "    #fig = plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.plot(sig)\n",
    "    #plt.scatter(locs,x[locs],marker = 'o',color = 'red')\n",
    "    #if len(c_locs) != 0:\n",
    "        #plt.scatter(c_locs,x[c_locs],marker = 'o',color = 'blue')\n",
    "\n",
    "    if len(c_locs) != 0:\n",
    "        ibis = c_locs\n",
    "    else:\n",
    "        ibis = locs\n",
    "\n",
    "    return ibis\n",
    "\n",
    "def d2max(x, fs):\n",
    "    \"\"\"\n",
    "    Detects inter-beat intervals using D2Max\n",
    "    Citation: Elgendi M, Norton I, Brearley M, Abbott D, Schuurmans D (2013) Systolic Peak Detection in Acceleration\n",
    "              Photoplethysmograms Measured from Emergency Responders in Tropical Conditions. PLoS ONE, vol. 8, no. 10,\n",
    "              pp. e76585. DOI: 10.1371/journal.pone.0076585\n",
    "\n",
    "    Inputs:   x, pulsatile signal [user defined units]\n",
    "              fs, sampling rate [Hz]\n",
    "    Outputs:  ibis, position of the starting points of inter-beat intervals [number of samples]\n",
    "\n",
    "    Developed by:  Elisa Mejía Mejía\n",
    "                   City, University of London\n",
    "    Version:       1.0 -   June, 2022\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Bandpass filter\n",
    "    if len(x) < 4098:\n",
    "        z_fill = np.zeros(4098 - len(x) + 1)\n",
    "        x_z = np.append(x, z_fill)\n",
    "    sos = sp.butter(10, [0.5, 8], btype = 'bp', analog = False, output = 'sos', fs = fs)\n",
    "    x_f = sp.sosfiltfilt(sos, x_z)\n",
    "\n",
    "    # Signal clipping\n",
    "    ind, = np.where(x_f < 0)\n",
    "    x_c = x_f\n",
    "    x_c[ind] = 0\n",
    "\n",
    "    # Signal squaring\n",
    "    x_s = x_c**2\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.plot(x_z)\n",
    "    #plt.plot(x_f)\n",
    "    #plt.plot(x_c)\n",
    "    #plt.plot(x_s)\n",
    "\n",
    "    # Blocks of interest\n",
    "    w1 = (111e-3)*fs\n",
    "    w1 = int(2*np.floor(w1/2) + 1)\n",
    "    b = (1/w1)*np.ones(w1)\n",
    "    ma_pk = sp.filtfilt(b,1,x_s)\n",
    "\n",
    "    w2 = (667e-3)*fs\n",
    "    w2 = int(2*np.floor(w2/2) + 1)\n",
    "    b = (1/w2)*np.ones(w1)\n",
    "    ma_bpm = sp.filtfilt(b,1,x_s)\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x_s/np.max(x_s))\n",
    "    #plt.plot(ma_pk/np.max(ma_pk))\n",
    "    #plt.plot(ma_bpm/np.max(ma_bpm))\n",
    "\n",
    "    # Thresholding\n",
    "    alpha = 0.02*np.mean(ma_pk)\n",
    "    th_1 = ma_bpm + alpha\n",
    "    th_2 = w1\n",
    "    boi = (ma_pk > th_1).astype(int)\n",
    "\n",
    "    blocks_init, = np.where(np.diff(boi) > 0)\n",
    "    blocks_init = blocks_init + 1\n",
    "    blocks_end, = np.where(np.diff(boi) < 0)\n",
    "    blocks_end = blocks_end + 1\n",
    "    if blocks_init[0] > blocks_end[0]:\n",
    "        blocks_init = np.append(1, blocks_init)\n",
    "    if blocks_init[-1] > blocks_end[-1]:\n",
    "        blocks_end = np.append(blocks_end, len(x_s))\n",
    "    #print('Initial locs BOI: ' + str(blocks_init))\n",
    "    #print('Final locs BOI: ' + str(blocks_end))\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x_s[range(len(x))]/np.max(x_s))\n",
    "    #plt.plot(boi[range(len(x))])\n",
    "\n",
    "    # Search for peaks inside BOIs\n",
    "    len_blks = np.zeros(len(blocks_init))\n",
    "    ibis = np.zeros(len(blocks_init))\n",
    "    for i in range(len(blocks_init)):\n",
    "        ind, = np.where(blocks_end > blocks_init[i])\n",
    "        ind = ind[0]\n",
    "        len_blks[i] = blocks_end[ind] - blocks_init[i]\n",
    "        if len_blks[i] >= th_2:\n",
    "            aux = x[blocks_init[i]:blocks_end[ind]]\n",
    "            if len(aux) != 0:\n",
    "                max_val = np.max(aux)\n",
    "                max_ind, = np.where(max_val == aux)\n",
    "                ibis[i] = max_ind + blocks_init[i] - 1\n",
    "\n",
    "    ind, = np.where(len_blks < th_2)\n",
    "    if len(ind) != 0:\n",
    "        for i in range(len(ind)):\n",
    "            boi[blocks_init[i]:blocks_end[i]] = 0\n",
    "    ind, = np.where(ibis == 0)\n",
    "    ibis = (np.delete(ibis, ind)).astype(int)\n",
    "\n",
    "    #plt.plot(boi[range(len(x))])\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.scatter(ibis, x[ibis], marker = 'o',color = 'red')\n",
    "\n",
    "    return ibis\n",
    "\n",
    "def upslopes(x):\n",
    "    \"\"\"\n",
    "    Detects inter-beat intervals using Upslopes\n",
    "    Citation: Arguello Prada EJ, Serna Maldonado RD (2018) A novel and low-complexity peak detection algorithm for\n",
    "              heart rate estimation from low-amplitude photoplethysmographic (PPG) signals. J Med Eng Technol, vol. 42,\n",
    "              no. 8, pp. 569-577. DOI: 10.1080/03091902.2019.1572237\n",
    "\n",
    "    Inputs:   x, pulsatile signal [user defined units]\n",
    "    Outputs:  ibis, position of the starting points of inter-beat intervals [number of samples]\n",
    "\n",
    "    Developed by:  Elisa Mejía Mejía\n",
    "                   City, University of London\n",
    "    Version:       1.0 -   June, 2022\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Peak detection\n",
    "    th = 6\n",
    "    pks = np.empty(0)\n",
    "    pos_pk = np.empty(0)\n",
    "    pos_pk_b = 0\n",
    "    n_pos_pk = 0\n",
    "    n_up = 0\n",
    "\n",
    "    for i in range(1, len(x)):\n",
    "        if x[i] > x[i - 1]:\n",
    "            n_up = n_up + 1\n",
    "        else:\n",
    "            if n_up > th:\n",
    "                pos_pk = np.append(pos_pk, i)\n",
    "                pos_pk_b = 1\n",
    "                n_pos_pk = n_pos_pk + 1\n",
    "                n_up_pre = n_up\n",
    "            else:\n",
    "                pos_pk = pos_pk.astype(int)\n",
    "                #print('Possible peaks: ' + str(pos_pk) + ', number of peaks: ' + str(n_pos_pk))\n",
    "                if pos_pk_b == 1:\n",
    "                    if x[i - 1] > x[pos_pk[n_pos_pk - 1]]:\n",
    "                        pos_pk[n_pos_pk - 1] = i - 1\n",
    "                    else:\n",
    "                        pks = np.append(pks, pos_pk[n_pos_pk - 1])\n",
    "                    th = 0.6*n_up_pre\n",
    "                    pos_pk_b = 0\n",
    "            n_up = 0\n",
    "    ibis = pks.astype(int)\n",
    "    #print(ibis)\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.scatter(ibis, x[ibis], marker = 'o',color = 'red')\n",
    "\n",
    "    return ibis\n",
    "\n",
    "def delineator(x, fs):\n",
    "    \"\"\"\n",
    "    Detects inter-beat intervals using Delineator\n",
    "    Citation: Li BN, Dong MC, Vai MI (2010) On an automatic delineator for arterial blood pressure waveforms. Biomed\n",
    "    Signal Process Control, vol. 5, no. 1, pp. 76-81. DOI: 10.1016/j.bspc.2009.06.002\n",
    "\n",
    "    Inputs:   x, pulsatile signal [user defined units]\n",
    "              fs, sampling rate [Hz]\n",
    "    Outputs:  ibis, position of the starting points of inter-beat intervals [number of samples]\n",
    "\n",
    "    Developed by:  Elisa Mejía Mejía\n",
    "                   City, University of London\n",
    "    Version:       1.0 -   June, 2022\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    # Lowpass filter\n",
    "    od = 3\n",
    "    sos = sp.butter(od, 25, btype = 'low', analog = False, output = 'sos', fs = fs)\n",
    "    x_f = sp.sosfiltfilt(sos, x)\n",
    "    x_m = 1000*x_f\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x)\n",
    "    #plt.plot(x_f)\n",
    "    #plt.plot(x_m)\n",
    "\n",
    "    # Moving average\n",
    "    n = 5\n",
    "    b = (1/n)*np.ones(n)\n",
    "    x_ma = sp.filtfilt(b,1,x_m)\n",
    "\n",
    "    # Compute differentials\n",
    "    dif = np.diff(x_ma)\n",
    "    dif = 100*np.append(dif[0], dif)\n",
    "    dif_ma = sp.filtfilt(b,1,dif)\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x_ma)\n",
    "    #plt.plot(dif_ma)\n",
    "\n",
    "    # Average thresholds in original signal\n",
    "    x_len = len(x)\n",
    "    if x_len > 12*fs:\n",
    "        n = 10\n",
    "    elif x_len > 7*fs:\n",
    "        n = 5\n",
    "    elif x_len > 4*fs:\n",
    "        n = 2\n",
    "    else:\n",
    "        n = 1\n",
    "    #print(n)\n",
    "\n",
    "    max_min = np.empty(0)\n",
    "    if n > 1:\n",
    "        #plt.figure()\n",
    "        #plt.plot(x_ma)\n",
    "        n_int = np.floor(x_len/(n + 2))\n",
    "        #print('Length of intervals: ' + str(n_int))\n",
    "        for j in range(n):\n",
    "            # Searches for max and min in 1 s intervals\n",
    "            amp_min, ind_min, amp_max, ind_max = seek_local(x_ma, int(j*n_int), int(j*n_int + fs))\n",
    "            #plt.scatter(ind_min, amp_min, marker = 'o', color = 'red')\n",
    "            #plt.scatter(ind_max, amp_max, marker = 'o', color = 'green')\n",
    "            max_min = np.append(max_min, (amp_max - amp_min))\n",
    "        max_min_avg = np.mean(max_min)\n",
    "        #print('Local max and min: ' + str(max_min) + ', average amplitude: ' + str(max_min_avg))\n",
    "    else:\n",
    "        amp_min, ind_min , amp_max, ind_max = seek_local(x_ma, int(close_win), int(x_len))\n",
    "        #plt.figure()\n",
    "        #plt.plot(x_ma)\n",
    "        #plt.scatter(ind_min, amp_min, marker = 'o', color = 'red')\n",
    "        #plt.scatter(ind_max, amp_max, marker = 'o', color = 'green')\n",
    "        max_min_avg = amp_max - amp_min\n",
    "        #print('Local max and min: ' + str(max_min) + ', average amplitude: ' + str(max_min_avg))\n",
    "\n",
    "    max_min_lt = 0.4*max_min_avg\n",
    "\n",
    "    # Seek pulse beats by min-max method\n",
    "    step_win = 2*fs       # Window length to look for peaks/onsets\n",
    "    close_win = np.floor(0.1*fs)\n",
    "                          # Value of what is considered too close\n",
    "\n",
    "    pks = np.empty(0)     # Location of peaks\n",
    "    ons = np.empty(0)     # Location of onsets\n",
    "    dic = np.empty(0)     # Location of dicrotic notches\n",
    "\n",
    "    pk_index = -1          # Number of peaks found\n",
    "    on_index = -1          # Number of onsets found\n",
    "    dn_index = -1          # Number of dicrotic notches found\n",
    "\n",
    "    i = int(close_win)    # Initializes counter\n",
    "    while i < x_len:      # Iterates through the signal\n",
    "        #print('i: ' + str(i))\n",
    "        amp_min = x_ma[i] # Gets the initial value for the minimum amplitude\n",
    "        amp_max = x_ma[i] # Gets the initial value for the maximum amplitude\n",
    "\n",
    "        ind = i           # Initializes the temporal location of the index\n",
    "        aux_pks = i       # Initializes the temporal location of the peak\n",
    "        aux_ons = i       # Initializes the temporal location of the onset\n",
    "\n",
    "        # Iterates while ind is lower than the length of the signal\n",
    "        while ind < x_len - 1:\n",
    "            #print('Ind: ' + str(ind))\n",
    "            # Verifies if no peak has been found in 2 seconds\n",
    "            if (ind - i) > step_win:\n",
    "                #print('Peak not found in 2 s')\n",
    "                ind = i   # Refreshes the temporal location of the index\n",
    "                max_min_avg = 0.6*max_min_avg  # Refreshes the threshold for the amplitude\n",
    "                # Verifies if the threshold is lower than the lower limit\n",
    "                if max_min_avg <= max_min_lt:\n",
    "                    max_min_avg = 2.5*max_min_lt # Refreshes the threshold\n",
    "                break\n",
    "\n",
    "            # Verifies if the location is a candidate peak\n",
    "            if (dif_ma[ind - 1]*dif_ma[ind + 1]) <= 0:\n",
    "                #print('There is a candidate peak')\n",
    "                # Determines initial and end points of a window to search for local peaks and onsets\n",
    "                if (ind + 5) < x_len:\n",
    "                    i_stop = ind + 5\n",
    "                else:\n",
    "                    i_stop = x_len - 1\n",
    "                if (ind - 5) >= 0:\n",
    "                    i_start = ind - 5\n",
    "                else:\n",
    "                    i_start = 0\n",
    "\n",
    "                # Checks for artifacts of saturated or signal loss\n",
    "                if (i_stop - ind) >= 5:\n",
    "                    for j in range(ind, i_stop):\n",
    "                        if dif_ma[j] != 0:\n",
    "                            break\n",
    "                    if j == i_stop:\n",
    "                        #print('Artifact')\n",
    "                        break\n",
    "\n",
    "                # Candidate onset\n",
    "                #print('Looking for candidate onsets...')\n",
    "                #plt.figure()\n",
    "                #plt.plot(x_ma)\n",
    "                if dif_ma[i_start] < 0:\n",
    "                    if dif_ma[i_stop] > 0:\n",
    "                        aux_min, ind_min, _, _ = seek_local(x_ma, int(i_start), int(i_stop))\n",
    "                        #plt.scatter(ind_min, aux_min, marker = 'o', color = 'red')\n",
    "                        if np.abs(ind_min - ind) <= 2:\n",
    "                            amp_min = aux_min\n",
    "                            aux_ons = ind_min\n",
    "                #print('Candidate onset: ' + str([ind_min, amp_min]))\n",
    "                # Candidate peak\n",
    "                #print('Looking for candidate peaks...')\n",
    "                if dif_ma[i_start] > 0:\n",
    "                    if dif_ma[i_stop] < 0:\n",
    "                        _, _, aux_max, ind_max = seek_local(x_ma, int(i_start), int(i_stop))\n",
    "                        #plt.scatter(ind_max, aux_max, marker = 'o', color = 'green')\n",
    "                        if np.abs(ind_max - ind) <= 2:\n",
    "                            amp_max = aux_max\n",
    "                            aux_pks = ind_max\n",
    "                #print('Candidate peak: ' + str([ind_max, amp_max]))\n",
    "                # Verifies if the amplitude of the pulse is larger than 0.4 times the mean value:\n",
    "                #print('Pulse amplitude: ' + str(amp_max - amp_min) + ', thresholds: ' +\n",
    "                #      str([0.4*max_min_avg, 2*max_min_avg]))\n",
    "                if (amp_max - amp_min) > 0.4*max_min_avg:\n",
    "                    #print('Expected amplitude of pulse')\n",
    "                    # Verifies if the amplitude of the pulse is lower than 2 times the mean value:\n",
    "                    if (amp_max - amp_min) < 2*max_min_avg:\n",
    "                        #print('Expected duration of pulse')\n",
    "                        if aux_pks > aux_ons:\n",
    "                            #print('Refining onsets...')\n",
    "                            # Refine onsets:\n",
    "                            aux_min = x_ma[aux_ons]\n",
    "                            temp_ons = aux_ons\n",
    "                            for j in range(aux_pks, aux_ons + 1, -1):\n",
    "                                if x_ma[j] < aux_min:\n",
    "                                    aux_min = x_ma[j]\n",
    "                                    temp_ons = j\n",
    "                            amp_min = aux_min\n",
    "                            aux_ons = temp_ons\n",
    "\n",
    "                            # If there is at least one peak found before:\n",
    "                            #print('Number of previous peaks: ' + str(pk_index + 1))\n",
    "                            if pk_index >= 0:\n",
    "                                #print('There were previous peaks')\n",
    "                                #print('Duration of ons to peak interval: ' + str(aux_ons - pks[pk_index]) +\n",
    "                                #     ', threshold: ' + str([3*close_win, step_win]))\n",
    "                                # If the duration of the pulse is too short:\n",
    "                                if (aux_ons - pks[pk_index]) < 3*close_win:\n",
    "                                    #print('Too short interbeat interval')\n",
    "                                    ind = i\n",
    "                                    max_min_avg = 2.5*max_min_lt\n",
    "                                    break\n",
    "                                # If the time difference between consecutive peaks is longer:\n",
    "                                if (aux_pks - pks[pk_index]) > step_win:\n",
    "                                    #print('Too long interbeat interval')\n",
    "                                    pk_index = pk_index - 1\n",
    "                                    on_index = on_index - 1\n",
    "                                    #if dn_index > 0:\n",
    "                                    #    dn_index = dn_index - 1\n",
    "                                # If there are still peaks, add the new peak:\n",
    "                                if pk_index >= 0:\n",
    "                                    #print('There are still previous peaks')\n",
    "                                    pk_index = pk_index + 1\n",
    "                                    on_index = on_index + 1\n",
    "                                    pks = np.append(pks, aux_pks)\n",
    "                                    ons = np.append(ons, aux_ons)\n",
    "                                    #print('Peaks: ' + str(pks))\n",
    "                                    #print('Onsets: ' + str(ons))\n",
    "\n",
    "                                    tf = ons[pk_index] - ons[pk_index - 1]\n",
    "\n",
    "                                    to = np.floor(fs/20)\n",
    "                                    tff = np.floor(0.1*tf)\n",
    "                                    if tff < to:\n",
    "                                        to = tff\n",
    "                                    to = pks[pk_index - 1] + to\n",
    "\n",
    "                                    te = np.floor(fs/20)\n",
    "                                    tff = np.floor(0.5*tf)\n",
    "                                    if tff < te:\n",
    "                                        te = tff\n",
    "                                    te = pks[pk_index - 1] + te\n",
    "\n",
    "                                    #tff = seek_dicrotic(dif_ma[to:te])\n",
    "                                    #if tff == 0:\n",
    "                                    #    tff = te - pks[pk_index - 1]\n",
    "                                    #    tff = np.floor(tff/3)\n",
    "                                    #dn_index = dn_index + 1\n",
    "                                    #dic[dn_index] = to + tff\n",
    "\n",
    "                                    ind = ind + close_win\n",
    "                                    break\n",
    "                            # If it is the first peak:\n",
    "                            if pk_index < 0:\n",
    "                                #print('There were no previous peaks')\n",
    "                                pk_index = pk_index + 1\n",
    "                                on_index = on_index + 1\n",
    "                                pks = np.append(pks, aux_pks)\n",
    "                                ons = np.append(ons, aux_ons)\n",
    "                                #print('Peaks: ' + str(pks))\n",
    "                                #print('Onsets: ' + str(ons))\n",
    "                                ind = ind + close_win\n",
    "                                break\n",
    "\n",
    "            ind = ind + 1\n",
    "        i = int(ind + 1)\n",
    "\n",
    "    if len(pks) == 0:\n",
    "        return -1\n",
    "    else:\n",
    "        x_len = len(pks)\n",
    "        temp_p = np.empty(0)\n",
    "        for i in range(x_len):\n",
    "            temp_p = np.append(temp_p, pks[i] - od)\n",
    "        ttk = temp_p[0]\n",
    "        if ttk < 0:\n",
    "            temp_p[0] = 0\n",
    "        pks = temp_p\n",
    "\n",
    "        x_len = len(ons)\n",
    "        temp_o = np.empty(0)\n",
    "        for i in range(x_len):\n",
    "            temp_o = np.append(temp_o, ons[i] - od)\n",
    "        ttk = temp_o[0]\n",
    "        if ttk < 0:\n",
    "            temp_o[0] = 0\n",
    "        ons = temp_o\n",
    "\n",
    "    pks = pks + 5\n",
    "    ibis = pks.astype(int)\n",
    "\n",
    "    return ibis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Fiducial Point functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following functions are required to detect fiducial points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fiducial_points(x,pks,fs,vis):\n",
    "    \"\"\"\n",
    "    Description: Pulse detection and correction from pulsatile signals\n",
    "    Inputs:  x, array with pulsatile signal [user defined units]\n",
    "             pks, array with the position of the peaks [number of samples]\n",
    "             fs, sampling rate of signal [Hz]\n",
    "             vis, visualisation option [True, False]\n",
    "    Outputs: fidp, dictionary with the positions of several fiducial points for the cardiac cycles [number of samples]\n",
    "\n",
    "    Fiducial points:  1: Systolic peak (pks)\n",
    "                      2: Onset, as the minimum before the systolic peak (ons)\n",
    "                      3: Onset, using the tangent intersection method (ti)\n",
    "                      4: Diastolic peak (dpk)\n",
    "                      5: Maximum slope (m1d)\n",
    "                      6: a point from second derivative PPG (a2d)\n",
    "                      7: b point from second derivative PPG (b2d)\n",
    "                      8: c point from second derivative PPG (c2d)\n",
    "                      9: d point from second derivative PPG (d2d)\n",
    "                      10: e point from second derivative PPG (e2d)\n",
    "                      11: p1 from the third derivative PPG (p1)\n",
    "                      12: p2 from the third derivative PPG (p2)\n",
    "\n",
    "    Libraries: NumPy (as np), SciPy (Signal, as sp), Matplotlib (PyPlot, as plt)\n",
    "\n",
    "    Version: 1.0 - June 2022\n",
    "\n",
    "    Developed by: Elisa Mejía-Mejía\n",
    "                   City, University of London\n",
    "\n",
    "    Edited by: Peter Charlton (see \"Added by PC\")\n",
    "\n",
    "    \"\"\"\n",
    "    # First, second and third derivatives\n",
    "    d1x = sp.savgol_filter(x, 9, 5, deriv = 1)\n",
    "    d2x = sp.savgol_filter(x, 9, 5, deriv = 2)\n",
    "    d3x = sp.savgol_filter(x, 9, 5, deriv = 3)\n",
    "\n",
    "    #plt.figure()\n",
    "    #plt.plot(x/np.max(x))\n",
    "    #plt.plot(d1x/np.max(d1x))\n",
    "    #plt.plot(d2x/np.max(d2x))\n",
    "    #plt.plot(d3x/np.max(d3x))\n",
    "\n",
    "    # Search in time series: Onsets between consecutive peaks\n",
    "    ons = np.empty(0)\n",
    "    for i in range(len(pks) - 1):\n",
    "        start = pks[i]\n",
    "        stop = pks[i + 1]\n",
    "        ibi = x[start:stop]\n",
    "        #plt.figure()\n",
    "        #plt.plot(ibi, color = 'black')\n",
    "        aux_ons, = np.where(ibi == np.min(ibi))\n",
    "        if len(aux_ons) > 1:\n",
    "            aux_ons = aux_ons[0]\n",
    "        ind_ons = aux_ons.astype(int)\n",
    "        ons = np.append(ons, ind_ons + start)\n",
    "        #plt.plot(ind_ons, ibi[ind_ons], marker = 'o', color = 'red')\n",
    "    ons = ons.astype(int)\n",
    "    #print('Onsets: ' + str(ons))\n",
    "    #plt.figure()\n",
    "    #plt.plot(x, color = 'black')\n",
    "    #plt.scatter(pks, x[pks], marker = 'o', color = 'red')\n",
    "    #plt.scatter(ons, x[ons], marker = 'o', color = 'blue')\n",
    "\n",
    "    # Search in time series: Diastolic peak and dicrotic notch between consecutive onsets\n",
    "    dia = np.empty(0)\n",
    "    dic = np.empty(0)\n",
    "    for i in range(len(ons) - 1):\n",
    "        start = ons[i]\n",
    "        stop = ons[i + 1]\n",
    "        ind_pks, = np.intersect1d(np.where(pks < stop), np.where(pks > start))\n",
    "        ind_pks = pks[ind_pks]\n",
    "        ibi_portion = x[ind_pks:stop]\n",
    "        ibi_2d_portion = d2x[ind_pks:stop]\n",
    "        #plt.figure()\n",
    "        #plt.plot(ibi_portion/np.max(ibi_portion))\n",
    "        #plt.plot(ibi_2d_portion/np.max(ibi_2d_portion))\n",
    "        aux_dic, _ = sp.find_peaks(ibi_2d_portion)\n",
    "        aux_dic = aux_dic.astype(int)\n",
    "        aux_dia, _ = sp.find_peaks(-ibi_2d_portion)\n",
    "        aux_dia = aux_dia.astype(int)\n",
    "        if len(aux_dic) != 0:\n",
    "            ind_max, = np.where(ibi_2d_portion[aux_dic] == np.max(ibi_2d_portion[aux_dic]))\n",
    "            aux_dic_max = aux_dic[ind_max]\n",
    "            if len(aux_dia) != 0:\n",
    "                nearest = aux_dia - aux_dic_max\n",
    "                aux_dic = aux_dic_max\n",
    "                dic = np.append(dic, (aux_dic + ind_pks).astype(int))\n",
    "                #plt.scatter(aux_dic, ibi_portion[aux_dic]/np.max(ibi_portion), marker = 'o')\n",
    "                ind_dia, = np.where(nearest > 0)\n",
    "                aux_dia = aux_dia[ind_dia]\n",
    "                nearest = nearest[ind_dia]\n",
    "                if len(nearest) != 0:\n",
    "                    ind_nearest, = np.where(nearest == np.min(nearest))\n",
    "                    aux_dia = aux_dia[ind_nearest]\n",
    "                    dia = np.append(dia, (aux_dia + ind_pks).astype(int))\n",
    "                    #plt.scatter(aux_dia, ibi_portion[aux_dia]/np.max(ibi_portion), marker = 'o')\n",
    "                    #break\n",
    "            else:\n",
    "                dic = np.append(dic, (aux_dic_max + ind_pks).astype(int))\n",
    "                #plt.scatter(aux_dia, ibi_portion[aux_dia]/np.max(ibi_portion), marker = 'o')\n",
    "    dia = dia.astype(int)\n",
    "    dic = dic.astype(int)\n",
    "    #plt.scatter(dia, x[dia], marker = 'o', color = 'orange')\n",
    "    #plt.scatter(dic, x[dic], marker = 'o', color = 'green')\n",
    "\n",
    "    # Search in D1: Maximum slope point\n",
    "    m1d = np.empty(0)\n",
    "    for i in range(len(ons) - 1):\n",
    "        start = ons[i]\n",
    "        stop = ons[i + 1]\n",
    "        ind_pks, = np.intersect1d(np.where(pks < stop), np.where(pks > start))\n",
    "        ind_pks = pks[ind_pks]\n",
    "        ibi_portion = x[start:ind_pks]\n",
    "        ibi_1d_portion = d1x[start:ind_pks]\n",
    "        #plt.figure()\n",
    "        #plt.plot(ibi_portion/np.max(ibi_portion))\n",
    "        #plt.plot(ibi_1d_portion/np.max(ibi_1d_portion))\n",
    "        aux_m1d, _ = sp.find_peaks(ibi_1d_portion)\n",
    "        aux_m1d = aux_m1d.astype(int)\n",
    "        if len(aux_m1d) != 0:\n",
    "            ind_max, = np.where(ibi_1d_portion[aux_m1d] == np.max(ibi_1d_portion[aux_m1d]))\n",
    "            aux_m1d_max = aux_m1d[ind_max]\n",
    "            if len(aux_m1d_max) > 1:\n",
    "                aux_m1d_max = aux_m1d_max[0]\n",
    "            m1d = np.append(m1d, (aux_m1d_max + start).astype(int))\n",
    "            #plt.scatter(aux_m1d, ibi_portion[aux_dic]/np.max(ibi_portion), marker = 'o')\n",
    "            #break\n",
    "    m1d = m1d.astype(int)\n",
    "    #plt.scatter(m1d, x[m1d], marker = 'o', color = 'purple')\n",
    "\n",
    "    # Search in time series: Tangent intersection points\n",
    "    tip = np.empty(0)\n",
    "    for i in range(len(ons) - 1):\n",
    "        start = ons[i]\n",
    "        stop = ons[i + 1]\n",
    "        ibi_portion = x[start:stop]\n",
    "        ibi_1d_portion = d1x[start:stop]\n",
    "        ind_m1d, = np.intersect1d(np.where(m1d < stop), np.where(m1d > start))\n",
    "        ind_m1d = m1d[ind_m1d] - start\n",
    "        #plt.figure()\n",
    "        #plt.plot(ibi_portion/np.max(ibi_portion))\n",
    "        #plt.plot(ibi_1d_portion/np.max(ibi_1d_portion))\n",
    "        #plt.scatter(ind_m1d, ibi_portion[ind_m1d]/np.max(ibi_portion), marker = 'o')\n",
    "        #plt.scatter(ind_m1d, ibi_1d_portion[ind_m1d]/np.max(ibi_1d_portion), marker = 'o')\n",
    "        aux_tip = np.round(((ibi_portion[0] - ibi_portion[ind_m1d])/ibi_1d_portion[ind_m1d]) + ind_m1d)\n",
    "        aux_tip = aux_tip.astype(int)\n",
    "        tip = np.append(tip, (aux_tip + start).astype(int))\n",
    "        #plt.scatter(aux_tip, ibi_portion[aux_tip]/np.max(ibi_portion), marker = 'o')\n",
    "        #break\n",
    "    tip = tip.astype(int)\n",
    "    #plt.scatter(tip, x[tip], marker = 'o', color = 'aqua')\n",
    "\n",
    "    # Search in D2: A, B, C, D and E points\n",
    "    a2d = np.empty(0)\n",
    "    b2d = np.empty(0)\n",
    "    c2d = np.empty(0)\n",
    "    d2d = np.empty(0)\n",
    "    e2d = np.empty(0)\n",
    "    for i in range(len(ons) - 1):\n",
    "        start = ons[i]\n",
    "        stop = ons[i + 1]\n",
    "        ibi_portion = x[start:stop]\n",
    "        ibi_1d_portion = d1x[start:stop]\n",
    "        ibi_2d_portion = d2x[start:stop]\n",
    "        ind_m1d = np.intersect1d(np.where(m1d > start),np.where(m1d < stop))\n",
    "        ind_m1d = m1d[ind_m1d]\n",
    "        #plt.figure()\n",
    "        #plt.plot(ibi_portion/np.max(ibi_portion))\n",
    "        #plt.plot(ibi_1d_portion/np.max(ibi_1d_portion))\n",
    "        #plt.plot(ibi_2d_portion/np.max(ibi_2d_portion))\n",
    "        aux_m2d_pks, _ = sp.find_peaks(ibi_2d_portion)\n",
    "        aux_m2d_ons, _ = sp.find_peaks(-ibi_2d_portion)\n",
    "        # a point:\n",
    "        ind_a, = np.where(ibi_2d_portion[aux_m2d_pks] == np.max(ibi_2d_portion[aux_m2d_pks]))\n",
    "        ind_a = aux_m2d_pks[ind_a]\n",
    "        if (ind_a < ind_m1d):\n",
    "            a2d = np.append(a2d, ind_a + start)\n",
    "            #plt.scatter(ind_a, ibi_2d_portion[ind_a]/np.max(ibi_2d_portion), marker = 'o')\n",
    "            # b point:\n",
    "            ind_b = np.where(ibi_2d_portion[aux_m2d_ons] == np.min(ibi_2d_portion[aux_m2d_ons]))\n",
    "            ind_b = aux_m2d_ons[ind_b]\n",
    "            if (ind_b > ind_a) and (ind_b < len(ibi_2d_portion)):\n",
    "                b2d = np.append(b2d, ind_b + start)\n",
    "                #plt.scatter(ind_b, ibi_2d_portion[ind_b]/np.max(ibi_2d_portion), marker = 'o')\n",
    "        # e point:\n",
    "        ind_e, = np.where(aux_m2d_pks > ind_m1d - start)\n",
    "        aux_m2d_pks = aux_m2d_pks[ind_e]\n",
    "        ind_e, = np.where(aux_m2d_pks < 0.6*len(ibi_2d_portion))\n",
    "        ind_e = aux_m2d_pks[ind_e]\n",
    "        if len(ind_e) >= 1:\n",
    "            if len(ind_e) >= 2:\n",
    "                ind_e = ind_e[1]\n",
    "            e2d = np.append(e2d, ind_e + start)\n",
    "            #plt.scatter(ind_e, ibi_2d_portion[ind_e]/np.max(ibi_2d_portion), marker = 'o')\n",
    "            # c point:\n",
    "            ind_c, = np.where(aux_m2d_pks < ind_e)\n",
    "            if len(ind_c) != 0:\n",
    "                ind_c_aux = aux_m2d_pks[ind_c]\n",
    "                ind_c, = np.where(ibi_2d_portion[ind_c_aux] == np.max(ibi_2d_portion[ind_c_aux]))\n",
    "                ind_c = ind_c_aux[ind_c]\n",
    "                if len(ind_c) != 0:\n",
    "                    c2d = np.append(c2d, ind_c + start)\n",
    "                    #plt.scatter(ind_c, ibi_2d_portion[ind_c]/np.max(ibi_2d_portion), marker = 'o')\n",
    "            else:\n",
    "                aux_m1d_ons, _ = sp.find_peaks(-ibi_1d_portion)\n",
    "                ind_c, = np.where(aux_m1d_ons < ind_e)\n",
    "                ind_c_aux = aux_m1d_ons[ind_c]\n",
    "                if len(ind_c) != 0:\n",
    "                    ind_c, = np.where(ind_c_aux > ind_b)\n",
    "                    ind_c = ind_c_aux[ind_c]\n",
    "                    if len(ind_c) > 1:\n",
    "                        ind_c = ind_c[0]\n",
    "                    c2d = np.append(c2d, ind_c + start)\n",
    "                    #plt.scatter(ind_c, ibi_2d_portion[ind_c]/np.max(ibi_2d_portion), marker = 'o')\n",
    "            # d point:\n",
    "            if len(ind_c) != 0:\n",
    "                ind_d = np.intersect1d(np.where(aux_m2d_ons < ind_e), np.where(aux_m2d_ons > ind_c))\n",
    "                if len(ind_d) != 0:\n",
    "                    ind_d_aux = aux_m2d_ons[ind_d]\n",
    "                    ind_d, = np.where(ibi_2d_portion[ind_d_aux] == np.min(ibi_2d_portion[ind_d_aux]))\n",
    "                    ind_d = ind_d_aux[ind_d]\n",
    "                    if len(ind_d) != 0:\n",
    "                        d2d = np.append(d2d, ind_d + start)\n",
    "                        #plt.scatter(ind_d, ibi_2d_portion[ind_d]/np.max(ibi_2d_portion), marker = 'o')\n",
    "                else:\n",
    "                    ind_d = ind_c\n",
    "                    d2d = np.append(d2d, ind_d + start)\n",
    "                    #plt.scatter(ind_d, ibi_2d_portion[ind_d]/np.max(ibi_2d_portion), marker = 'o')\n",
    "    a2d = a2d.astype(int)\n",
    "    b2d = b2d.astype(int)\n",
    "    c2d = c2d.astype(int)\n",
    "    d2d = d2d.astype(int)\n",
    "    e2d = e2d.astype(int)\n",
    "    #plt.figure()\n",
    "    #plt.plot(d2x, color = 'black')\n",
    "    #plt.scatter(a2d, d2x[a2d], marker = 'o', color = 'red')\n",
    "    #plt.scatter(b2d, d2x[b2d], marker = 'o', color = 'blue')\n",
    "    #plt.scatter(c2d, d2x[c2d], marker = 'o', color = 'green')\n",
    "    #plt.scatter(d2d, d2x[d2d], marker = 'o', color = 'orange')\n",
    "    #plt.scatter(e2d, d2x[e2d], marker = 'o', color = 'purple')\n",
    "\n",
    "    # Search in D3: P1 and P2 points\n",
    "    p1p = np.empty(0)\n",
    "    p2p = np.empty(0)\n",
    "    for i in range(len(ons) - 1):\n",
    "        start = ons[i]\n",
    "        stop = ons[i + 1]\n",
    "        ibi_portion = x[start:stop]\n",
    "        ibi_1d_portion = d1x[start:stop]\n",
    "        ibi_2d_portion = d2x[start:stop]\n",
    "        ibi_3d_portion = d3x[start:stop]\n",
    "        ind_b = np.intersect1d(np.where(b2d > start),np.where(b2d < stop))\n",
    "        ind_b = b2d[ind_b]\n",
    "        ind_c = np.intersect1d(np.where(c2d > start),np.where(c2d < stop))\n",
    "        ind_c = c2d[ind_c]\n",
    "        ind_d = np.intersect1d(np.where(d2d > start),np.where(d2d < stop))\n",
    "        ind_d = d2d[ind_d]\n",
    "        ind_dic = np.intersect1d(np.where(dic > start),np.where(dic < stop))\n",
    "        ind_dic = dic[ind_dic]\n",
    "        #plt.figure()\n",
    "        #plt.plot(ibi_portion/np.max(ibi_portion))\n",
    "        #plt.plot(ibi_1d_portion/np.max(ibi_1d_portion))\n",
    "        #plt.plot(ibi_2d_portion/np.max(ibi_2d_portion))\n",
    "        #plt.plot(ibi_3d_portion/np.max(ibi_3d_portion))\n",
    "        #plt.scatter(ind_b - start, ibi_3d_portion[ind_b - start]/np.max(ibi_3d_portion), marker = 'o')\n",
    "        #plt.scatter(ind_c - start, ibi_3d_portion[ind_c - start]/np.max(ibi_3d_portion), marker = 'o')\n",
    "        #plt.scatter(ind_d - start, ibi_3d_portion[ind_d - start]/np.max(ibi_3d_portion), marker = 'o')\n",
    "        #plt.scatter(ind_dic - start, ibi_3d_portion[ind_dic - start]/np.max(ibi_3d_portion), marker = 'o')\n",
    "        aux_p3d_pks, _ = sp.find_peaks(ibi_3d_portion)\n",
    "        aux_p3d_ons, _ = sp.find_peaks(-ibi_3d_portion)\n",
    "        # P1:\n",
    "        if (len(aux_p3d_pks) != 0 and len(ind_b) != 0):\n",
    "            ind_p1, = np.where(aux_p3d_pks > ind_b - start)\n",
    "            if len(ind_p1) != 0:\n",
    "                ind_p1 = aux_p3d_pks[ind_p1[0]]\n",
    "                p1p = np.append(p1p, ind_p1 + start)\n",
    "                #plt.scatter(ind_p1, ibi_3d_portion[ind_p1]/np.max(ibi_3d_portion), marker = 'o')\n",
    "        # P2:\n",
    "        if (len(aux_p3d_ons) != 0 and len(ind_c) != 0 and len(ind_d) != 0):\n",
    "            if ind_c == ind_d:\n",
    "                ind_p2, = np.where(aux_p3d_ons > ind_d - start)\n",
    "                ind_p2 = aux_p3d_ons[ind_p2[0]]\n",
    "            else:\n",
    "                ind_p2, = np.where(aux_p3d_ons < ind_d - start)\n",
    "                ind_p2 = aux_p3d_ons[ind_p2[-1]]\n",
    "            if len(ind_dic) != 0:\n",
    "                aux_x_pks, _ = sp.find_peaks(ibi_portion)\n",
    "                if ind_p2 > ind_dic - start:\n",
    "                    ind_between = np.intersect1d(np.where(aux_x_pks < ind_p2), np.where(aux_x_pks > ind_dic - start))\n",
    "                else:\n",
    "                    ind_between = np.intersect1d(np.where(aux_x_pks > ind_p2), np.where(aux_x_pks < ind_dic - start))\n",
    "                if len(ind_between) != 0:\n",
    "                    ind_p2 = aux_x_pks[ind_between[0]]\n",
    "            p2p = np.append(p2p, ind_p2 + start)\n",
    "            #plt.scatter(ind_p2, ibi_3d_portion[ind_p2]/np.max(ibi_3d_portion), marker = 'o')\n",
    "    p1p = p1p.astype(int)\n",
    "    p2p = p2p.astype(int)\n",
    "    #plt.figure()\n",
    "    #plt.plot(d3x, color = 'black')\n",
    "    #plt.scatter(p1p, d3x[p1p], marker = 'o', color = 'green')\n",
    "    #plt.scatter(p2p, d3x[p2p], marker = 'o', color = 'orange')\n",
    "\n",
    "    # Added by PC: Magnitudes of second derivative points\n",
    "    bmag2d = np.zeros(len(b2d))\n",
    "    cmag2d = np.zeros(len(b2d))\n",
    "    dmag2d = np.zeros(len(b2d))\n",
    "    emag2d = np.zeros(len(b2d))\n",
    "    for beat_no in range(0,len(d2d)):\n",
    "        bmag2d[beat_no] = d2x[b2d[beat_no]]/d2x[a2d[beat_no]]\n",
    "        cmag2d[beat_no] = d2x[c2d[beat_no]]/d2x[a2d[beat_no]]\n",
    "        dmag2d[beat_no] = d2x[d2d[beat_no]]/d2x[a2d[beat_no]]\n",
    "        emag2d[beat_no] = d2x[e2d[beat_no]]/d2x[a2d[beat_no]]\n",
    "\n",
    "     # Added by PC: Refine the list of fiducial points to only include those corresponding to beats for which a full set of points is available\n",
    "    off = ons[1:]\n",
    "    ons = ons[:-1]\n",
    "    if pks[0] < ons[0]:\n",
    "        pks = pks[1:]\n",
    "    if pks[-1] > off[-1]:\n",
    "        pks = pks[:-1]\n",
    "\n",
    "    # Visualise results\n",
    "    if vis == True:\n",
    "        fig, (ax1,ax2,ax3,ax4) = plt.subplots(4, 1, sharex = True, sharey = False, figsize=(10,10))\n",
    "        fig.suptitle('Fiducial points')\n",
    "\n",
    "        ax1.plot(x, color = 'black')\n",
    "        ax1.scatter(pks, x[pks.astype(int)], color = 'orange', label = 'pks')\n",
    "        ax1.scatter(ons, x[ons.astype(int)], color = 'green', label = 'ons')\n",
    "        ax1.scatter(off, x[off.astype(int)], marker = '*', color = 'green', label = 'off')\n",
    "        ax1.scatter(dia, x[dia.astype(int)], color = 'yellow', label = 'dia')\n",
    "        ax1.scatter(dic, x[dic.astype(int)], color = 'blue', label = 'dic')\n",
    "        ax1.scatter(tip, x[tip.astype(int)], color = 'purple', label = 'dic')\n",
    "        ax1.legend()\n",
    "        ax1.set_ylabel('x')\n",
    "\n",
    "        ax2.plot(d1x, color = 'black')\n",
    "        ax2.scatter(m1d, d1x[m1d.astype(int)], color = 'orange', label = 'm1d')\n",
    "        ax2.legend()\n",
    "        ax2.set_ylabel('d1x')\n",
    "\n",
    "        ax3.plot(d2x, color = 'black')\n",
    "        ax3.scatter(a2d, d2x[a2d.astype(int)], color = 'orange', label = 'a')\n",
    "        ax3.scatter(b2d, d2x[b2d.astype(int)], color = 'green', label = 'b')\n",
    "        ax3.scatter(c2d, d2x[c2d.astype(int)], color = 'yellow', label = 'c')\n",
    "        ax3.scatter(d2d, d2x[d2d.astype(int)], color = 'blue', label = 'd')\n",
    "        ax3.scatter(e2d, d2x[e2d.astype(int)], color = 'purple', label = 'e')\n",
    "        ax3.legend()\n",
    "        ax3.set_ylabel('d2x')\n",
    "\n",
    "        ax4.plot(d3x, color = 'black')\n",
    "        ax4.scatter(p1p, d3x[p1p.astype(int)], color = 'orange', label = 'p1')\n",
    "        ax4.scatter(p2p, d3x[p2p.astype(int)], color = 'green', label = 'p2')\n",
    "        ax4.legend()\n",
    "        ax4.set_ylabel('d3x')\n",
    "\n",
    "        plt.subplots_adjust(left = 0.1,\n",
    "                            bottom = 0.1,\n",
    "                            right = 0.9,\n",
    "                            top = 0.9,\n",
    "                            wspace = 0.4,\n",
    "                            hspace = 0.4)\n",
    "\n",
    "    # Creation of dictionary\n",
    "    fidp = {'pks': pks.astype(int),\n",
    "            'ons': ons.astype(int),\n",
    "            'off': off.astype(int),  # Added by PC\n",
    "            'tip': tip.astype(int),\n",
    "            'dia': dia.astype(int),\n",
    "            'dic': dic.astype(int),\n",
    "            'm1d': m1d.astype(int),\n",
    "            'a2d': a2d.astype(int),\n",
    "            'b2d': b2d.astype(int),\n",
    "            'c2d': c2d.astype(int),\n",
    "            'd2d': d2d.astype(int),\n",
    "            'e2d': e2d.astype(int),\n",
    "            'bmag2d': bmag2d,\n",
    "            'cmag2d': cmag2d,\n",
    "            'dmag2d': dmag2d,\n",
    "            'emag2d': emag2d,\n",
    "            'p1p': p1p.astype(int),\n",
    "            'p2p': p2p.astype(int)\n",
    "            }\n",
    "\n",
    "    return fidp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
